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Abstract 


This  report  summarizes  methodologies  and  findings  of  a  study  performed  at  the  School  of  Engineering, 
University  of  Glasgow,  funded  by  the  European  Office  of  Aerospace  Research  and  Development,  during  the  period 
March  2013  to  August  2014. 

The  study  aimed  at  developing  algorithms  for  accurate,  long  term  propagation  of  a  particular  type  of  space  debris, 
made  of  multi-layer  insulation  (MLI)  foils,  originated  from  delamination  from  ageing  spacecraft.  This  type  of  debris 
has  high  area-to-mass  ratio,  combined  with  high  reflectivity  index,  and  hence  it  is  subject  to  strong  solar  radiation 
pressure  (SRP)  acceleration.  In  addition,  MLI  membranes  are  highly  flexible,  therefore  their  effective,  exposed  area 
to  the  sun  is  subject  to  change  over  time.  This  effect  has  been  taken  into  account  in  the  past  through  averaging 
techniques.  This  study  attempts,  for  the  first  time,  to  include  the  dynamics  of  the  deformation  in  the  propagation  of 
the  equations  of  motion  (orbit  and  attitude).  The  computational  cost  of  modelling  a  deformable  membrane  subject  to 
forces  can  be  significant,  hence  two  different  models  were  developed:  one  uses  a  linear  Bernoulli -Euler  beam  model, 
while  the  second  uses  discretizes  the  membrane  properties  through  lump  masses,  spring  and  damper.  The  latter 
framework  models  arbitrarily  large  deformations,  as  expected  on  the  low-bending-stiffness  membrane. 

Propagations  of  a  typical  geostationary  orbit  (GEO)  debris  using  both  models  are  compared  to  other  models  not 
taking  into  account  the  deformation  (cannonball  model  or  rigid  membrane).  The  presented  results  show  that 
considerable  difference  in  the  estimation  of  the  orbital  parameters  (particularly  inclination  and  eccentricity)  can  be 
obtained  even  over  tens  of  days. 

In  addition,  in  order  to  reduce  the  computational  burden  of  the  equations  of  motion,  a  numerical  integrator  has 
been  developed,  being  able  to  treat  coupled  dynamics  with  different  typical  time -scales,  specifically  orbit  and 
attitude  equations.  Assuming  the  attitude  is  the  fast  dynamics,  it  resulted  that  faster  integration  times  can  be  obtained, 
with  the  same  level  of  accuracy  of  a  traditional  algorithm,  if  the  computational  cost  of  the  fast  dynamics  is  above  a 
certain  threshold. 

An  appendix  shows  the  validation  of  the  deformation  computed  through  the  linear  model  with  the  commercial 
FEM  solver  ABAQUS. 
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1.  Introduction 


Debris  objects  are  uncontrolled,  non-functional  objects  orbiting  around  the  Earth.  Natural  debris  includes 
asteroids  and  rocks,  while  artificial,  human-made  debris  comprise  non-active  spacecraft  (including  launcher  upper 
stages)  and  their  components,  fragmented  following  collision  events  or  separation  due  to  ageing.  Debris  objects  often 
travel  at  very  high  velocities,  and  for  this  reason,  despite  their  size  can  be  relatively  small,  they  that  can  threat  active 
spacecraft,  leading  to  catastrophic  break-ups,  which  in  turn  generate  new  space  debris.  If  the  rate  of  production  from 
collision  is  greater  than  decay  rate  (for  example  due  to  atmospheric  reentry),  then  a  chain  reaction  starts,  which  leads 
to  the  growth  of  the  number  of  debris  objects  exponentially,  as  predicted  by  Kessler.  This  scenario  is  now  known  as 
“Kessler  syndrome”. 

Since  the  early  space  era,  the  launch  of  Sputnik  on  4  October  1957,  more  than  5,000  satellites  and  upper  stages 
were  launched  between  low  Earth  orbit  (LEO)  and  geostationary  Earth  orbit  (GEO)  and  led  to  increasing  population 
of  more  than  22,000  trackable  objects  with  sizes  larger  than  10  cm.  Approximately  1,000  of  these  are  operational 
spacecraft.  The  remaining  94%  are  space  debris,  i.e  defunct  satellites  and  rocket  upper  stages  which  are  not  useful 
purpose.  The  routinely  tracked  objects  (about  64%)  are  fragments  from  some  250  breakups,  mainly  explosions  and 
collisions  of  satellites  or  rocket  bodies  and  space  experiment.  An  estimated  700,000  objects  larger  than  1  cm  and  170 
million  objects  larger  than  1  mm  are  expected  to  be  in  the  Earth  orbits.  Fig.  1  (credit:  NASA)  shows  a  snapshot  of 
known  debris  objects  (not  to  scale)  in  LEO  and  GEO  range. 


Fig.  1.  Objects  in  Earth  orbit  that  are  currently  being  tracked.  Approximately  95%  of  the  objects  in  this 
illustration  are  orbital  debris  [Souce:  NASA]. 


One  of  the  most  important  debris -causing  events  in  history  happened  on  11  January  2007.  A  direct-ascent, 
kinetic-kill  anti-satellite  (AS AT)  vehicle  was  successfully  tested  by  destroying  an  inactive  Chinese  Feng  Yun  1C 
(FY-1C)  weather  satellite  [1].  This  event  alone  generated  more  debris  objects  than  any  previous  space  accident.  The 
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US  Space  Surveillance  Network  (SSN)  was  able  to  detect  more  than  900  pieces  of  debris  that  were  at  risk  to  damage 
operational  spacecraft. 

In  February  10,  2009,  defunct  Russian  military  satellite  Cosmos  2251  crashed  the  active  U.S.  communications 
satellite  Iridium  33  at  790  km  above  Siberia  [2].  This  collision  created  two  large  debris  clouds  and  the  SSN  reported 
that  382  pieces  of  debris  from  Iridium  33  and  893  pieces  from  Cosmos  2251  were  created,  and  they  will  possibly 
remain  in  orbit  for  up  to  10,000  years. 

Fig.  2  shows  some  evidence  of  damage  made  by  space  debris  to  active  spacecraft. 


b) 

Fig.  2.  (a)  Hole  made  by  debris  in  the  panel  of  NASA’s  Solar  Max  experiment,  (b)  Window  of  the  space  shuttle 
damaged  by  a  tiny  piece  of  space  junk  (a  paint  fleck)  during  the  STS -7  mission  (Photo:  NASA  Orbital  Debris 
Program  Office) 

1.1.  Survey  on  orbital  debris  and  models 

1.1.1.  Fragmentation  and  distribution  models 

Based  on  history  of  collision  and  explosions  of  spacecraft  and  upper  stages  and  the  observation  data  in  low  Earth 
orbit,  the  NASA  Orbital  Debris  Program  Office,  established  since  the  1970s,  studied  and  implemented  a  debris 
fragmentation  distribution  model  which  was  implemented  into  EVOLVE.  The  model  implements  an  estimation  of 
size,  area-to-mass  and  velocity  of  debris  and  the  current  version  is  EVOLVE  4.0. 

Anz-Meador  [3]  validated  debris  sizes  in  the  range  1  mm  to  10  cm,  which  pose  major  hazard  to  spacecraft.  Both 
explosion  models  from  EVOLVE  and  LEGEND  were  utilized  to  observe  on-orbit  fragmentation  debris  arising  from 
the  explosion  of  rocket  bodies  and  payloads  and  ground-based  laboratory  tests.  The  data  products  for  validation  are 
from  the  US  SSN,  NASA  Haystack  and  Auxiliary  (HAX)  radar  observation  campaigns  and  Satellite  Orbital  Debris 
Characterization  Impact  Test  (SOCIT).  It  resulted  that  the  number  and  size  of  the  distribution  is  against  the  observed 
environment  and  conclusions  were  drawn  that  further  tests  were  required  before  adopting  the  model.  Instead,  the  area 
to  mass  distribution  is  adequate  to  describe  the  ballistic  properties  of  both  large  and  small  debris  objects. 
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When  the  number  of  fragments  resulting  from  a  collision  is  of  the  order  of  hundreds  or  thousands,  numerical 
simulations  of  their  propagation  become  computationally  very  expensive,  up  to  a  point  that  they  are  not  viable.  For 
this  reason,  Li  Yi-Yong  [4]  studied  and  developed  novel  approaches  to  propagate  debris  clouds  accurately. 

EVOLVE  was  originally  capable  to  simulate  the  orbit  debris  in  LEO,  however  it  was  expanded  into 
GEO_EVOLVE  1.0  to  simulate  the  orbital  debris  in  the  Earth  geosynchronous  region  since  1999.  The  models  were 
to  be  updated  because  of  the  overestimation  of  the  future  collision  rate  in  GEO,  due  to  variation  of  argument  of 
perigee  and  precession  of  the  nodes.  The  next  update  of  GEO_EVOLVE  2.0  [5]  improved  two  significant  points: 
firstly,  the  development  traffic  database  was  upgraded  with  the  inclusion  of  reliable  insertion  node,  argument  of 
perigee  and  mean  anomaly  of  all  launches.  Secondly,  perturbations  were  added,  such  as  atmospheric  drag,  lunar 
gravitational  perturbation,  solar  radiation  pressure,  Earth  gravity  field  zonal  (J2,  J3  and  J4)  and  tesseral  (J2,2,  J34,  etc.) 
harmonics.  GEOPROP  and  PROP3D  propagators  were  used. 

Another  LEO-to-GEO  evolutionary  debris  model,  LEGEND  [6],  was  developed  by  the  NASA  Debris  Program 
Office  at  Johnson  Space  Center.  It  is  based  on  historical  satellite  launch  database  (DBS -database)  to  reproduce  the 
debris  environment  between  1957  and  2001  and  predict  future  debris  environment.  It  covers  the  near  Earth  space 
between  200  and  2,000  km,  medium  Earth  orbit  (MEO)  and  GEO  (34,000-38,000  km)  and  super  GEO  (above  38,000 
km).  This  model  provided  all  debris  characteristics.  The  two  propagators  used  in  LEGEND  were  GEOPROP  and 
PROP3D.  The  GEOPROP  propagator  calculated  the  motion  of  GEO  objects  while  the  PROP3D  propagator  focused 
on  the  orbits  of  all  other  satellites  in  the  simulation. 

All  models  are  investigated  from  NASA,  While  the  European  Space  Agency  (ESA)  started  to  develop  the 
Meteoriod  And  Space  Debris  Terrestrial  Environment  Reference  model  (MASTER)  in  1995  [7].  The  model  was 
introduced  to  derive  spatial  densities  and  velocity  distributions  in  a  3D  control  volume  from  LEO  and  GEO  altitude 
for  the  debris  size  larger  than  1  mm,  similarly  to  NASA’s  EVOLVE  and  LEGEND.  MASTER  2009  added  data  from 
new  collision  and  observations  by  the  TIRA,  EISCAT,  Goldstone  and  Haystack  radars  and  the  ESA  Space  Debris 
Telescope  (ESA-SDT).  High  area-to-mass  ratio  (HAMR)  objects  were  also  considered  [8].  According  to  the 
observed  data  by  ESA  space  debris  telescope  to  characterize  the  space  debris  environment  near  the  geostationary 
ring  [9],  there  were  2,790  uncorrelated  targets  brighter  than  magnitude  18.5.  A  Gaussian  distribution  was  applied  to 
approximate  the  number  of  detections,  focusing  on  objects  with  magnitude  higher  than  18.5  with  a  probability  about 
95%.  The  study  determined  detection  probabilities  for  a  given  field  of  view  (FOV)  and  compared  the  ESA  data:  the 
results  were  in  good  agreement  with  the  debris  population  models  of  MASTER. 

1.1.2.  Debris  in  GTO 

In  geosynchronous  transfer  orbits  (GTO),  debris  are  mainly  created  by  disposed  rocket  upper  stages  and 
fragments  from  rocket  explosions  and  collisions  [10].  Most  of  GTOs  perigees  are  in  LEO,  below  2000  km,  and 
apogees  near  or  above  GEO.  Therefore,  debris  in  GTO  has  the  potential  for  interfering  with  objects  in  both  LEO  and 
GEO,  which  is  the  reason  why  strategies  and  approaches  were  suggested  to  minimize  space  debris  in  GTO  [11],  like 
re-entry  at  perigee  in  LEO.  In  a  study  of  the  probability  of  collision  (PC)  of  spacecraft  in  GEO,  Mcknight  [12] 


11 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


0!^  afcs.« ®  B0  University  i  School  of 

Dynamics  of  flexible  MLI-type  debris  for  accurate  orbit  prediction  '  'J*MJ  m  -  -S>«^  -  of  Glasgow  Engineering 

considered  parameters  as  spatial  density  (SPD),  relative  collision  velocity  (VR),  collision  cross-section  of  satellite  at 
risk  (AC)  and  time  (T). 

1.1.3.  High  area-to-mass  ratio  debris 

The  observation  by  Schildknecht  [13]  with  the  ESA  1 -metre  telescope  (ESASDT)  [14]  and  Zimmerwald  Laser 
and  Astrometry  Telescope  (ZIMLAT)  revealed  small-size  debris  objects  in  highly  elliptical  orbits.  It  found  the  new 
type  of  debris  population  for  7  unknown  which  was  small  size  of  debris  in  elliptical  orbits.  Three  of  the  identified 
objects  orbited  in  a  classical  GTO  with  eccentricity  about  0.7.  Other  four,  instead,  had  semi-major  axes  (and  hence 
period)  close  to  that  of  the  GEO  [15],  but  the  range  of  eccentricities  was  between  0.13  and  0.49.  Subsequently, 
around  80  objects  of  the  new  population  were  found  but  the  original  sources  were  unknown,  until  it  was  assumed 
that  they  were  pieces  of  multi-layer  insulation  (MLI)  material.  This  material  is  very  lightweight  (high  area-to-mass 
ratio)  and  highly  reflective:  the  solar  radiation  pressure  significantly  perturbs  HAMR  objects,  resulting  in  periodical 
variations  of  eccentricity  and  inclination.  In  that  work,  the  numerical  integrator  SATORB  was  applied  to  with  a  force 
model  including  the  geopotential  harmonics  up  order  twelve.  Other  perturbations  included  the  SRP,  the  solar  and 
lunar  gravitational  forces  and  the  Earth’s  shadow  effects.  Lastly,  the  objects  were  propagated  assuming  random 
tumbling  motion  and  shapes.  The  light  curves  showed  a  wide  variety  of  signatures,  ranging  from  periodic  or  random 
variations  of  several  magnitudes  over  time  spans  of  a  few  minutes  to  constant  brightness  over  10  min.  It  was  also 
observed  degradation  of  FEP  material  properties  of  the  MLI  around  the  telescope  Hubble  space  telescope  and  this 
material  was  returned  to  Earth  for  deep  investigation  [16].  Today,  there  is  suspect  that  this  type  of  debris  is  made  of 
multi-layer  insulation,  but  there  is  no  evidence  to  confirm  so. 

Further  experiments  were  conducted  by  Murakami  et  al.  to  investigate  fragmentation  of  multi-layer  insulation 
[17],  and  results  were  compared  with  NASA  breakup  model,  showing  similar  trend  in  term  of  mass  distribution. 
Some  evidence  supports  that  the  mass  distribution  concerning  material  in  the  experiments  of  the  study  of  Murakami 
[17]  and  support  that  most  of  High- Area-mass-ratio  were  MLI  that  was  one  of  debris  in  space.  Anselmo  et  al.  [18] 
investigated  the  long-term  evolution  of  HAMR  debris,  showing  variation  of  inclination  and  eccentricity  over 
centuries. 

A  novel  approach  to  determine  the  average  cross-sectional  area  of  breakup  fragments  was  proposed  by  Hanada  et 
al.  [19].  This  approach  provides  a  good  estimate  in  the  calculation  of  the  cross-sectional  area  of  multi-layer 
insulation  fragments. 

1.1.4 .  Effect  of  SRP 

Light  curves  are  often  used  to  determine  the  properties  of  debris  objects.  There  are  however  unresolved  light 
curves  of  objects  orbiting  in  GEO,  which  could  be  matched  with  MLI-type  of  objects.  The  investigation  in  [20] 
defined  probable  initial  conditions  of  angular  momentum,  torque,  quaternion  representation  of  orientation.  Linear 
combinations  of  the  database  were  considered  to  model  objects  with  both  specular  and  diffuse  reflection. 
Furthermore,  angular  momentum  and  torque  were  calculated  using  solar  radiation  pressure  considering  specular  and 
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diffuse  reflection.  The  three  types  of  MLI  were  investigated,  and  according  to  the  study,  the  angular  momentum 
depended  on  the  area-to-mass  ratio,  the  irregularity  of  shape  reflectivity  and  the  internal  friction. 

Flegel  [21]  investigated  MLI  and  rocket  upper  stages  and  parts  of  satellites  (LMRO)  in  term  of  total  spatial  object 
density  on  altitude  and  over  time  from  1958  to  2010.  The  number  of  MLI  objects  peaked  around  900-1400  km  and 
decreased  slowly  in  GEO  while  LMRO  population  distribution  is  peak  in  range  between  1,000  and  5,000.  In  the 
GEO  region,  the  peak  of  number  LMRO  was  at  25,000  km  and  36,000  km.  The  number  of  MLI-type  objects  was 
higher  than  that  of  LMRO.  The  model  was  validated  comparing  the  observation  data  from  ESA’s  PROOF  (Program 
for  Radar  and  Optical  Observation  Forcasting). 

Liou  et  al.  [23]  proposed  that  the  HARM  could  be  a  plausible  explanation  for  the  debris  with  highly  eccentric 
orbits  with  period  similar  to  GEO.  The  global  behavior  of  thousands  of  HAMR  debris  objects  revealed  the  unique 
pattern  of  their  distribution  in  space. 

The  orbit  determination  of  high  area-to-mass  ratio  objects  in  GEO  trajectories  was  analyzed  by  Kelecy  et  al.  [22]: 
the  cooler  cover  of  a  MSG-2  satellite  was  used  for  the  investigation  and  the  effects  of  the  non-modeled  accelerations 
was  considered  to  predict  the  HAMR  object  trajectory.  The  analyzed  objects  were  assumed  to  be  flat  plates,  covered 
with  MLI.  The  simulator  included  Earth  gravity  zonal  harmonics,  luni-solar  gravitational  perturbations  and  attitude 
dependent  solar  radiation  pressure  acceleration  that  considered  reflective,  diffusive  and  absorption  components.  An 
extended  Kalman  filter  (EKF)  was  used  to  estimate  the  orbit  in  the  orbit  determination  tool  kit  (ODTK).  The  initial 
position  and  velocity  states  were  derived  using  a  combination  of  Gauss -Gooding  initial  orbit  determination  (IOD). 

Kelecy  et  al.  [23]  studied  the  solar  radiation  pressure  effect  on  a  variety  of  high  area-to-mass  ratio  debris  in  GEO 
orbit  by  focusing  on  MLI  or  thermal  insulation,  to  improve  the  orbit  predictions.  It  was  found  that  the  characteristics 
and  dynamics  of  the  objects  could  be  identified  by  tracking  and  measuring  data. 

A  semi-analytical  theory  was  used  to  analyze  the  gravitational  perturbation  of  debris  in  GEO  by  Valk  et  al.  [24], 
who  also  considered  the  Earth’s  shadowing  effects  on  HAMR  objects  in  the  geosynchronous  region.  Because  of  the 
non-singular  set  of  equations  used,  the  method  was  particularly  suitable  for  both  near  circular  and  near  equatorial 
orbits  as  well  as  orbits  that  transit  periodically  around  null  eccentricities  and  null  inclinations.  The  short,  mid  and 
long  term  propagation  under  SRP  and  the  shadowing  effects  was  applied  to  objects  of  20  m2/kg  and  the  numerical 
and  semi-analytical  methods  were  compared  over  different  time  scales.  Inclination,  eccentricity  and  semi-major  axis 
were  used  to  assess  the  results  of  the  propagations. 

Another  investigation  of  the  shadow  effects  on  high  area-to-mass  geostationary  space  debris,  in  the  short  and 
long  term  time  evolution,  was  done  by  Hubaux  [25].  The  mean  exponential  growth  factor  of  nearby  orbits  (MEGNO) 
was  introduced  to  investigate  both  shadow  conical  and  cylindrical  Earth’s  shadowing  model.  It  resulted  that  the 
cylindrical  model  was  unreliable  for  predicting  the  orbital  behavior,  while  one  of  conical  modelled  was  proven  to  be 
better  because  of  its  realistic  umbra-penumbra  transition  model. 

MEGNO  was  applied  by  Valk  [26]  to  propagate  the  dynamics  of  HAMR  in  GEO  in  the  long  term.  The  validation 
of  the  method  was  by  propagation  on  over  12,600  orbits  and  over  30  years  of  time  span.  Perturbations  were 
considered  such  as  the  influence  of  the  Earth’s  gravity  field,  both  the  gravitational  perturbation  of  the  sun  and  the 
moon,  and  the  direct  solar  radiation  pressure.  A  sensitivity  analysis  to  initial  conditions  was  done.  Secondary 
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resonances  on  both  sides  of  the  well-known  pendulum  like  the  pattern  of  geostationary  objects  were  found.  In  this 
work,  MEGNO  proved  to  be  an  efficient  method  to  propagate  the  phase  space. 

Another  study  in  long-term  dynamics  of  HAMR  debris  is  in  [18]:  a  range  of  objects,  including  and  multi-layer 
insulation  (MLI)  blankets  and  other  low  density  composite  material,  were  considered.  The  selected  initial  conditions 
were  that  of  objects  released  in  geosynchronous  orbit,  and  a  high-accuracy  numerical  integrator  was  used,  which 
included  the  zonal  and  tesseral  harmonics  of  the  Earth,  shadow  effects  and  atmospheric  drag.  Interestingly,  the 
dimensionless  radiation  pressure  coefficient  Cr  was  considered  for  taking  into  account  the  deformed  geometry  of 
debris,  affecting  the  area  to  mass  ratio,  such  that  Cr  x  Aim  gives  the  effective  area-to-mass  ratio  exposed  to  the  SRP. 

The  resonance  effect  caused  by  the  geometry  of  the  Earth -Moon- Sun  systems  in  case  of  HAMR  objects  was  also 
studied  by  Rosengren  [27].  It  was  noted  that  the  mean  geometry  of  the  Earth-Moon-Sun  system  repeats  itself  closely 
after  a  period  of  time,  of  about  223  synodic  months,  resulting  in  a  dynamical  system  moving  in  a  nearly  periodic 
orbit.  When  displayed  in  the  heliocentric  orbit  frame,  the  actual  motion  resulted  very  complicate. 

In  many  studies,  pieces  of  multi-layer  insulation  blankets  are  suspected  to  be  HARM  debris  but  there  is  no 
observational  proof  of  this.  Most  “classical”  low  AMR  GEO  debris  objects  have  flat  or  slowly  varying  light  curves, 
instead  HAMR  objects  show  all  variations  with  period  of  a  few  minutes  and  shorter.  One  of  the  issues  is  that  the 
knowledge  of  the  surface  properties  and  the  shapes  of  the  objects  including  specular  and  diffuse  reflection,  albedo, 
and  color  are  required  to  identify  these  types  of  debris. 

The  observation  of  HAMR  objects  from  ESA  Space  Debris  Telescope  (ESASDT)  was  used  for  the  study  by  Fruh 
et  al.  [28],  who  investigated  the  variation  of  their  area-to-mass  ratio  and  orbital  parameters  and  showed  that 
shadowing  effects  were  significant  influence  on  the  long  term  evolution  of  the  orbits.  However,  in  this  work,  the 
AMR  values  were  assumed  to  be  constant:  the  results  were  not  the  exact  dynamics  of  HAMR  objects  because  effects 
due  to  irregular  shape,  complex  attitude  motion  and  deformation  of  the  actual  objects  were  not  considered.  These 
effects  produce  an  actual  change  of  the  AMR  value  over  time. 

More  recently,  a  further  investigation  of  MLI  debris  in  near  geosynchronous  orbits  was  done  by  Fruh,  [29].  Here, 
MLI  objects  were  assumed  to  be  flat  plates  with  material  properties  of  PET  and  Kapton,  as  these  materials  are  those 
mostly  used  to  protect  spacecraft.  Propagation  considered  gravitational  field,  direct  solar  radiation  pressure,  Earth 
magnetic  field  and  the  influence  of  shadow  effect  with  a  simple  cylindrical  shadow  model.  An  important  point  of  this 
study  is  that  it  coupled  both  orbit  and  attitude  propagation  over  4  days.  Then,  light  curve  were  used  to  correlate  the 
attitude.  The  results  demonstrated  that  the  solar  radiation  torque  and  gravitational  torque  led  to  rapid  attitude  motion, 
while  magnetic  torques  stabilize  the  attitude  motion  preventing  a  fast  spin.  Resulting  overestimation  of  the  evolution 
of  eccentricity  and  inclinations  originated  from  the  assumption  of  a  constant  AMR  value. 

Another  study  of  the  effect  of  SRP  on  orbital  debris  is  that  of  Fruh  and  Schildknecht  [28].  The  work  used  several 
objects,  including  an  upper  stage  rocket  debris,  Block  D-2,  simplified  in  a  cylindrical  shape,  and  other  HAMR 
objects  modelled  as  flat  rigid  sheets,  which  were  tracked  over  a  number  of  years.  The  attitude  motion  was  correlated 
with  the  observed  light  curves,  which  measured  the  variation  of  the  brightness  of  the  object  due  to  lighting  condition, 
viewing  geometry,  shape,  reflection  properties  and  attitude  motion.  The  considered  perturbations  were  the  third  body 
gravitational  field  and  solar  radiation  pressure.  The  orbit  attitude  motion  was  calculated  by  Euler’s  equation  using  the 
assumption  of  debris  as  rigid  body.  The  study  showed  that  the  flat  sheet  was  subject  to  fast  oscillations,  while  the 
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attitude  motion  of  the  upper  stage  was  significantly  influenced  by  an  offset  between  geometrical  center  and  the 
center  of  mass  object. 

The  study  of  the  evolution  of  the  relative  velocity  and  the  eccentricity  of  particles  due  to  SRP  was  done  by 
Celestino  et  al.  [10].  Different  particle  sizes  were  considered:  according  to  the  results,  particles  size  less  than  5  pm 
were  significantly  affected  by  SRP.  These  sizes  of  debris  experienced  a  very  large  amplitude  of  oscillation  in  the 
eccentricity  value  over  time,  leading  to  a  natural  mechanism  of  removal  when  passing  close  the  region  of  influence 
of  air  drag.  Instead,  the  SRP  affected  the  larger  debris  particles  in  the  short  period  inducing  changes  in  eccentricity 
and  semi-major  axis. 
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2.  Dynamics 


2.1.  Orbital  dynamics 

2.1.1.  Orbital  elements 

Orbital  elements  are  the  parameters  to  identify  a  specific  orbit  around  the  Earth.  The  common  sets  usually  are 
classical  orbital  elements,  or  Keplerian  elements  illustrated  in  Fig.  3. 

Z 


Fig.  3.  Orbital  elements. 


2. 1.1.1  Dimensional  Elements 

These  describe  the  shape  and  size  of  the  orbit  and  relate  the  position  on  orbit: 

•  A  is  the  semi-major  axis  of  the  orbit,  identifying  its  size. 

•  E  is  the  eccentricity  of  the  orbit  determining  how  elliptical  the  orbit  is. 

•  ns  the  time  of  the  satellite  measured  from  the  periapsis  passage  (closest  approach  to  the  main  attractor). 
It  relates  to  the  mean  motion,  M.  It  determines  the  satellite’s  current  position.  Alternatively,  the  true 
anomaly  0  can  be  used. 

2.1. 1.2  Orientation  Elements 

•  i  is  the  inclination,  i.e.  the  angle  between  the  orbit  plane  and  respect  to  the  reference  plane.  If  0°  <  i  < 
90°  ,  the  motion  is  termed  as  direct  or  prograde  and  retrograde  if  90°  <  i  <180°. 

•  co  is  the  argument  of  perigee,  i.e.  is  the  angle  from  the  ascending  node  to  the  eccentricity  vector, 
measured  anticlockwise. 

•  Q  is  the  right  ascension  of  the  ascending  node,  or  RAAN,  the  angle  from  the  direction  of  the  vernal 
equinox  to  the  line  of  the  nodes,  measured  counterclockwise. 


16 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


0!^  afcs.« ®  B0  University  i  School  of 

Dynamics  of  flexible  MLI-type  debris  for  accurate  orbit  prediction  '  d  -  -S>«^  -  of  Glasgow  Engineering 

The  six  classical  orbital  elements  are  useful  and  intuitive  to  describe  the  orbit  and  current  position  of  satellites. 
However,  the  equations  to  describe  the  dynamics  are  highly  nonlinear. 

2.1.2.  Coordinate  systems 

2. 1.2.1  Earth-centered  inertial  (ECI)  coordinate  system 

The  origin  is  located  in  the  center  of  the  Earth,  and  the  three  axes  point  into  inertial  directions.  The  z  axis  points 
towards  the  North  Pole  and  it  is  perpendicular  to  the  equatorial  plane.  The  x  axis  points  towards  the  vernal  equinox 
(first  point  of  Aries),  and  the  y  axis  completes  the  right  handed  frame. 

2. 1.2.1  Earth-centered  Earth-fixed  (ECEF)  coordinate  system 

The  origin  is  located  in  the  center  of  the  Earth,  and  the  three  axes  are  fixed  to  the  Earth.  The  z  axis  points  towards 
the  North  Pole  and  it  is  perpendicular  to  the  equatorial  plane.  The  x  axis  points  towards  the  Greenwich  meridian,  and 
the  y  axis  completes  the  right  handed  frame.  This  frame  is  used  to  describe  the  position  of  terrestrial.  Fig.  4  shows 
the  both  coordinate  systems. 


Fig.  4.  ECI  and  ECEF  coordinate  systems. 

2.1.3.  Two-body  problem 

Newton’s  law  of  universal  gravitation  states  that  every  point  mass  attracts  every  single  point  mass  by  a  force 
pointing  towards  the  other  mass  and  proportional  to  the  product  of  the  two  masses  and  inversely  proportional  to  the 
square  of  the  distance  between  them.  In  this  research,  we  focus  on  a  body  (debris), of  mass  m  ,  orbiting  around  the 
Earth,  whose  mass  is  m0  : 

mm  _ 


where  G  is  the  universal  gravitation  constant,  6.674  x  10'11  Nm2kg'2  ,  and  r  is  the  position  vector  of  the  object  from 
the  origin  in  the  center  of  the  Earth. 
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K 


Fig.  5.  Gravitational  force  from  the  Earth  acting  on  a  satellite 


By  introducing  the  approximation  that  the  orbiting  mass  is  much  smaller  than  that  of  the  Earth,  it  is  possible  to 
obtain  the  following  differential  equation  of  motion  for  the  object  with  respect  to  the  center  of  the  Earth: 

r  =  -  — r 


where  ju  =  G  (  ra  ffi  +  m  )  «  G m 


2.1.4 .  Third  body  perturbations 


For  an  object  orbiting  the  Earth,  the  effects  of  gravitational  attraction  from  other  (third)  bodies,  such  as  the  Sun 
and  Moon  become  more  significant  with  the  orbit  altitude.  In  geosynchronous  orbit,  the  main  third  body  gravitational 
perturbing  forces  acting  on  the  objects  are  the  Sun  and  Moon.  The  equation  of  motion,  including  the  gravitational 
perturbation  of  N  other  bodies  (as  derived  in  [30])  is: 


r 


(  r  )  N 

—  l+G2> 

Vr  ) 


where  rk  is  the  position  vector  of  the  k-th  attractor  (Sun,  Moon)  and  mk  its  mass. 


2.1.5.  Solar  radiation  pressure 

The  solar  radiation  pressure  (SRP)  is  caused  by  the  incident  rays  of  light  (photons)  from  the  Sun  on  the  exposed 
surfaces  of  the  object.  It  is  a  non-conservative  and  non-gravitational  disturbance.  This  perturbation  in  GEO  is  usually 
dominant,  especially  for  objects  which  have  large  surfaces  exposed  to  the  sun,  and  small  masses  (i.e.  high  area-to- 
mass  ratio).  Modelling  the  SRP  perturbing  acceleration  is  difficult  for  a  number  of  reasons. 

First  of  all,  the  solar  flux  generally  varies  depending  on  the  time  of  the  year  and  the  solar  activity,  which  has  its 
own  cycle.  In  addition,  for  an  accurate  model  of  the  acceleration,  all  the  areas  exposed  to  the  sun  have  to  be 
considered.  This  includes  modelling  the  Earth  eclipses  during  the  orbit.  Finally,  the  reflective  properties  of  the 
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materials  on  the  surfaces  exposed  to  the  sun  are  also  important  to  determine  the  force  applied  to  each  one  of  them. 

The  development  of  an  accurate  model  concerns  the  effectiveness  to  calculate  the  position  of  the  sun  with  respect  to 

the  object  and  the  Earth  at  any  given  time  during  the  orbit  as  well  as  the  correct  attitude  of  the  objects. 

According  to  Wertz  [31],  solar  flux  value  is  measured  at  Earth  aphelion  point  is: 

1358  /  2 

SF  =  - W /m 

1.004  +  0.03 34  cos  D  ,  , 


where  d  aphelion  is  defined  as  In  times  the  time  elapsed  from  Earth  aphelion,  expressed  as  a  fraction  of  the  whole  year. 

The  amount  of  momentum  that  is  carried  by  this  flux  can  be  found  through  Einstein’s  law,  e  =  me2 ,  and  can  be 
expressed  as  a  force  per  unit  area,  or  pressure,  as 

sf  6  N 

p  =  - =  4.57x10  - 

1  sr  2 

c  m 


Where  c  is  the  speed  of  light.  The  numerical  value  is  approximate  at  1  Astronomical  Unit  (AU). 

The  actual  value  and  direction  of  the  force  experienced  by  a  flat  surface  subject  to  psr  depends  on  the  orientation 
of  the  surface  with  respect  to  the  sun  and  the  reflectivity  properties  of  the  surface.  Fig.  6  shows  the  effect  of  specular 
reflection  of  the  photons  and  diffuse  reflection,  typical  of  a  real,  rough  surface.  In  addition,  Lambertian  diffusion  is 
considered  the  diffuse  and  specular  radiation  forces  that  affect  to  the  object  and  consider  each  object  surface  area  as 
following.  In  the  ideal  case  of  pure  specular  reflection,  the  force  applied  is  perpendicular  to  the  surface.  In  non-ideal 
cases,  the  force  is  deflected  by  a  certain  angle,  as  diffuse  reflection  and  Lambertian  diffusion  generate  tangential 
forces  (see  Fig.  7  for  reflection  geometry).  The  same  figure  shows  the  incident  angle,  defined  as  the  angle  between 
the  normal  unit  vector  to  the  surface  and  the  unit  vector  defining  the  direction  from  the  Sun  to  the  surface.  The 
reflected  angle  is  equal  to  the  incident  in  the  case  of  purely  specular  reflection. 

The  acceleration  due  to  SRP  can  be  computed,  considering  the  summation  over  each  flat  surface  i  exposed  to  the 
sun,  as: 


a 


rad 


"  P„A ,COS0  [  ( CRd  \  1 

Z  - "  2|  - ~+CRi  cos(f  )  |  n  +  (l  -  C  )s  j- 

IT,  m  U  3  '  j  J 


(1) 


where  A.  is  the  area  of  the  i-th  surface,  m  is  the  mass  of  the  entire  object,  and  c rd  and  c rs  are  the  coefficients  of 
diffuse  and  specular  reflectivity,  respectively  [30].  If  we  consider  that  the  surface  can  also  absorb  a  fraction  of  the 
photons,  given  by  a  coefficient  c  ra  ,  then  it  must  be  Crs  +  cds  #  cra  =  l  .  The  surface  normal  unit  vector,  n  ,  and  the 

sun  vector  s  are  required  to  be  known,  from  the  orientation  of  the  specify  the  orientation  of  the  object.  The 
summation  adds  the  N  flat  plate  surfaces  of  the  satellite  model. 
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Specular  Reflection  Diffuse  Reflection 


Fig.  6.  Specular  and  diffuse  reflection. 


Reflected  beam 

Fig.  7.  Incident  solar  radiation  geometry.  Note  that  the  unit  vectors  t  and  n  are  the  unit  vectors  in  the 
tangential  and  normal  directions  (with  respect  to  the  surface),  respectively;  h  and  in  are  oriented  upwards 
from  the  plane  of  the  paper;  s  and  rsato  are  collinear  with  the  incident  solar  rays,  pointing  towards  the  sun. 


2.1.6.  Shadow  effect 

The  shadow  effect  is  a  partial  or  total  obscuration  of  the  sunlight  from  the  satellite  by  the  Earth.  A  simplified 
cylindrical  Earth  shadow  model,  as  that  in  Fig.  8,  does  not  consider  penumbra:  the  object  is  either  in  sunlight  or 
complete  shadow.  A  conical  model,  instead,  uses  a  more  realistic  conical  Earth  shadow  with  two  main  regions, 
penumbra  and  umbra  as  in  Fig.  9.  In  the  penumbra  region,  the  sun  is  partially  eclipsed  from  the  Earth.  Clearly,  the 
solar  radiation  pressure  is  zero  completely  when  satellite  is  in  the  Earth  shadow.  It  is  instead  more  complex  to  find 
the  SRP  value  in  the  penumbra  region,  and  this  limits  the  applicability  of  this  model. 


satellite 


Exit 


✓  Earth 


!  Enter 


Fig.  8.  Cylindrical  eclipse  geometry. 
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Fig.  9.  Eclipse  geometry  for  conical  model. 


In  this  work,  we  will  the  cylindrical  model.  Whether  the  spacecraft  is  in  umbra  or  sunlight  can  be  computed 
according  to  Vallado  [30]  through  the  following  parameters: 


P-r-jh  !'+('■  •r2)i 


where  Re  =  6, 37 1  km  is  the  mean  radius  of  the  Earth.  Then, 

•  If  Tmin  <  o  v  zmin  >  o  v  |c  Tmin  f  >  l  ,  the  object  is  in  sunlight; 

•  Else,  the  object  is  in  umbra. 

The  SRP  acceleration  is  simply  set  to  zero  when  the  object  is  in  umbra,  otherwise  the  actual  value  is  used,  as 
computed  in  the  previous  section. 


2.1.7.  Sun  vector 

The  unit  vector  s  defines  the  direction  from  satellite  to  the  Sun  at  any  given  time  as: 

s  = 

where  rsatQ  is  the  vector  from  the  satellite  to  the  Sun  as  in  Fig.  10,  as: 

rsatO  =  rffi®  -  r 

According  to  Vallado  [30],  the  Earth-Sun  vector  reQ  is  calculated  at  any  given  time  by  using  the  following 

formulas.  First,  we  define  the  number  of  Julian  centuries  as  a  function  of  the  Julian  Date  JD: 

JD  -  2451545 
T  =  - 

TDB 

36525 

Then  the  mean  longitude  A M  e  ,  the  solar  mean  anomaly  m  0  ,  the  ecliptic  longitude  A ecllptic ,  the  distance  of  the  sun 
with  respect  to  the  Earth  rQ  ,  and  the  obliquity  of  the  ecliptic,  s,  from  the  following  equations: 
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M  Q  =  357.5277233°  +  35999. 05034T 


AmQ  =  28  0.460°  +  3  6000.77  0  TTj 


Aecl  =  AmQ  +  1.9  14666471  sin  (Mq  )  +  0.0 1 9994643  sin  (2  M  Q) 


rQ  =1.000140612  -  0 .0 1 6708 6  cos  ( M  G  )  -  0 .000 1 3 95  89  cos  (2 M  0  ) 
£  =  2  3.4  3  9  291° -  0.0  13  0042TrflB  -1.64xl0"7  T*DB  +  5 .04  x  1  0  "7  T*DB 

Given  kecliptic  and  s  at  epoch,  the  sun  position  vector  with  respect  to  the  Earth  is: 

r  r0  cos(Aecllpllc)  1 
rG  =  j  rG  cosO)  sin (kecliptic)  j 
[r0  sin(£)sin(Xec(ip<(c)  J 


where  all  distances  are  in  AU. 


Fig.  10.  Earth-Satellite-Sun  vector  diagram 


2.1.8.  Moon  vector 

The  Moon  position  with  respect  to  the  Earth  is  computed  in  a  similar  way  [30]: 

k  .  =  218.32°  +  481,  267. 883T  +  6.29  sin  (134.9  +  477, 198. 85T  )  - 

ecliptic  7  TDB  \  7  TDB  ) 

1.27  sin  (25  9.2  -  413,335.38 TJDB  )  +  0.66  sin  (235  .7  +  890, 5  34.23rrDB  )  + 
0.2  1  sin  (2  69.9  +  9  5  4, 3  97 .70TrDB  )  -  0.19  sin  (3  5  7.5  +  3  5,999.05 Ttdb  )  - 
0.1 1  sin  (18  6.6  +  9  6  6,404.05Trz)B  ) 

</>ecliptic  =  5.1  3°  (93.3  +  4  81,202.03rrDB  )  +  0.2  3  sin  (228.2  +  9  60, ,  400.8  7rrDB  )  - 
0.28  sin  (3  18.3  +  600  3.1  8TrDB  )  -  0 . 1  7  sin  ( 2  1 7 .6  +  407, 332.20Ttdb  ) 

p  =  0.95  08°  +  0.0518  sin  (134.9  +  477,19  8.8  5TrDB  )  + 

0.009  5  cos  (25  9.2  +  4  1  3, 3  35 .1  8TrDB  )  +  0.007  8  cos  (2  3  5.7  +  8  9  0,5  34.2  3TrDB  )  + 
0.002  8  sin  (269.9  +  95  4,  3 97 .7 0TrDB  ) 

1 

^ Moon 

s  i  n  p 
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where  Aeclipm ,  <l>ecliptic  are  the  ecliptic  longitude  and  latitude  of  the  Moon  respectively,  p  is  the  horizontal  parallax,  and 
rMoon  *s  distance  of  the  Moon  from  the  Earth.  Finally, 


r 


Moon 


1  cos  {d>  .  )  cos(/l  ...  )  \ 

|  ecliptic  7  v  ecliptic  '  | 

I  cos(f )  cos( (ft ec[iptic )  c o  s ( A, eciiptic )  —  s in ( £ )  sin ( $ eciiptic )  | 
sin(^ )cos(^  )cos(/l  +  cos(f)sin(^  ) 

|_  v  7  v“  ecliptic  7  v  ecliptic'  v  7  ecliptic  '  J 


2.2.  Attitude  dynamics 

The  angular  momentum  h  of  a  rigid  body  with  respect  to  its  center  of  mass  is  defined  as  by  the  following 
equation: 

h=f  rx(rxco)  dm 


r  denotes  the  location  vector  of  a  particle  mass  element  dm  and  to  is  the  angular  velocity  vector  of  body. 

The  total  torque  on  the  system  about  its  center  of  mass  equals  the  time  derivative  of  the  angular  momentum  in 
inertial  coordinates: 


N 


Or,  in  a  body-fixed  frame,  we  get  the  classical  Euler  attitude  equations: 


dh 

T  t  =  <j  -  (0  X  (  1(0  ) 


where 


body 


are  the  angular  momentum  vector  change  of  a  rigid  body  with  respect  to  the  center  of 


mass  in  an  inertial  frame  and  in  a  body  frame,  t  f  is  the  total  external  torque  acting  on  the  body,  I  is  the  matrix  of 
moment  of  inertia  of  the  rigid  body  and  «  is  the  angular  velocity  vector  in  the  body  frame. 

Note  that  these  equations  are  found  under  the  assumption  that  I  is  constant,  i.e.  a  rigid  body.  We  therefore  imply 
that  the  deformations  are  slow  with  respect  to  any  change  in  attitude,  so  changes  in  I  can  be  neglected. 

In  this  work,  the  modelled  torques  are:  the  solar  radiation  pressure  torque  t  rad  ,  the  gravity  gradient  t  Grav  and  the 

magnetic  torque  t  Mag  .  The  total  torque  is  then: 


T 


+  T 


2.2.7.  SRP  torque 

From  the  solar  radiation  pressure  force  Frad  computed  before,  the  torque  is  found  as: 

t  = p  . xTF 

rad  r  rad  rad 
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where  p  rad  is  the  vector  from  the  center  of  mass  of  the  object  to  the  center  of  pressure,  and  T  is  the  transformation 
matrix  from  the  geocentric  reference  frame  to  the  body  reference  system. 


2.2.2.  Gravity  gradient  torque 

The  gravitational  torque  for  a  rigid  body  is  derived  by  determining  the  gravitational  force  F  this  equation  can  be 
written  as  [32]: 

3 gm 

tg  =  - ^[TrxITr] 


2.2.3.  Magnetic  torque 

The  motion  of  a  conductor  (i.e.  a  metallic  object)  in  a  magnetic  field  (i.e.  that  of  the  Earth)  induces  Eddy 
currents,  creating  forces  and  torques  [31].  The  induced  torque  can  be  given  as: 

TMag  =  k  I®  xTB]x  TB 

where  B  is  the  local  Earth  magnetic  field,  is  approximated  in  polar  coordinates  as: 

f  R®  V  r  0  1  1  n 

B e  =  ^ - j  |  g1  sin#  -  gj  cos#cos^  -  hl  cos#  sin 

(  Rm  V  r  ,  !  n 

- J  [g1  sin  (j)  -  hx  cos^J 

(  R®  3  r  0  1  1  -I 

B  =  2  ^ - j  cos#  -  g  j  sin#cos^  -  hx  sin#sin^J 

where  the  coefficients  for  the  Earth  are:  g°  =  -30109  nT  ,  g\  =  -2006  nT  ,  h\  =  5704.2  nT  ,  and  6  and  ^  are  geodetic 
latitude  and  longitude  respectively.  Finally, 

n  3 

k  =  — a  l  A 
4 

with  /  the  radius  of  the  loop  on  the  area  of  surface  A,  o  the  conductivity  of  the  material. 


2.2.4.  Moments  of  inertia 

In  this  research,  as  a  preliminary  assumption  which  we  aim  to  remove  in  future  work,  and  only  for  the  purpose  of 
easing  the  calculation  of  the  matrix  of  inertia,  we  assume  that  the  debris  is  flat.  In  principal  axes  of  inertia,  the  matrix 
of  inertia  I  becomes: 

\IXX  0  0  1 

1  =  !  0  /„  0  j 

L°  0  uJ 

For  a  flat  plate  of  dimensions  w,  h  and  d : 
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ix 


i 

yy 


I 

zz 


=  — m  ( h 2  +  d  2  ) 
12  V 


=  — m  ( w  2  +  d  2 ) 
12 


=  — m  (  w 2  +  h 2  ) 
12 


2.3.  Rotational  kinematics 


The  attitude  of  a  tumbling  body  is  most  conveniently  integrated  with  a  set  of  axes  fixed  to  the  body  (principal 
axes  of  inertia).  It  is  therefore  required  to  parameterize  the  attitude  to  obtain  the  transformation  matrix  R  from  the 
body  axes  (superscript  B)  to  a  set  of  inertial  coordinates  (N): 


r  ru  r12 
\  R  21  R  22 
L^31  R32 


*13  1 
R23  \  yN 

R»\ 


The  rotation  matrix  R  can  be  expressed  as  a  function  of  pitch  (<p),  yaw  (y/)  and  roll  ( 6 )  angles,  by  using  a  3-1-3 
Euler  rotation,  as  [32]. 


f  cos  \|/  sin  \j/ 

D  BIN  I  . 

R  =  |  sin  \| /  cos  \| / 

|  0  0 


oin 

°!!° 

i  J  L° 


0  0  ~|  [cos  (p  sin  (p  0 ~| 

cos  6  sin  6  jj  sin  cp  cos  cp  0  j 
sin  6  cos$J[o  0  lj 


leading  to: 


r  c(\| f)c{(p)-  j(\|/)c(<9)5(^) 

Rb,n  =j  -c(\\f)s((p)- s(y)c(0)s((p) 

[  s(0) 


c(\\l)c(0)s{(p)+  s(\\l)c((p)  s(0)s(\ |/)1 

c(\\f)c(0)c{(p)-  5(\|/)^(^)  s{e)s{(p)  | 

-c(\i/)s(6»)  c(6)  \ 


where  cos()  is  abbreviated  by  “c()”  and  sin()  by  cts()”. 

While  Euler  angles  provide  an  intuitive  way  of  defining  the  attitude  at  any  given  time,  they  have  singularities  and 
therefore  they  are  not  suitable  for  being  used  in  automatic  integration. 


2.3.1.  Quaternions 

Quaternions  are  an  alternative  way  of  parameterizing  the  attitude,  with  no  singularities,  at  the  cost  of  using  four 
scalar  variables  instead  of  three.  It  is  possible  to  rotate  from  an  initial  orientation  to  final  orientation  with  a  single 
rigid  rotation,  around  an  eigenvector  e,  defining  a  principal  rotation  axis,  of  a  principal  rotation  angle  ®,  through  the 
matrix: 

R  =  cos  (<J>  )  1  +  |T  -  cos  (O  )J  [et  e2  ^3][^!  e3  ]  _sin(0)E 

r  0  -e3  e2  1 

with:  E  =  j  e3  0  -e1  j 

|_-e2  et  0  J 


25 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Dynamics  of  flexible  MLI-type  debris  for  accurate  orbit  prediction 


|g  University  School  of 
of Glasgow  Engineering 


or: 

6  xc  3  (l  -  cos  (®  ))  -  e2  sin  (®  )1 
e2e3(l-cos(®))  +  ^  sin  (®  )  | 
cos(®  )  +  e3  (1  -  cos(®  ))  j 

The  principal  axis  can  be  found  starting  from  R  as: 


2  sin  ® 

2  sin  ® 

2  sin  ® 


[  cos(0  )  +  e~  (1  -  cos(®  ))  exe2  (l  -  cos  (®  ))  +  e3  sin  (®  ) 
R  =  |  gjg2  (l- cos(0  ))-  e3  s  in  (®)  cos(®)  +  e2(l-cos(®)) 
j^g1e3(l-cos(®))+g2sin(®)  e2e3(l-cos(®))-^1sin(®) 


Finally,  the  quaternions  can  be  defined  as: 


® 

qx  =  e  sin( — ) 
2 

® 

q2  =  e2  sin( — ) 
2 

® 

q3  =  e3  sin( — ) 


a 

q4  =  cos(-) 
2 


or,  in  vector  notation:  q  = 


i  i 
_  i  1 


I  q  I 
I  I 
L^J 


The  singularity  is  removed  because  only  three  of  the  quaternions  are  actually  necessary  to  reconstruct  the 
attitude,  and  the  following  relationship,  instead  of  one  of  the  above,  can  be  used  to  find  the  value  of  the  singular 
quaternion: 

q  i  +  q  2  +  4  3  +  q]  =  1 


Finally,  the  cosine  matrix  to  transform  an  arbitrary  vector  can  be  expressed  directly  as  a  function  of  the 
quaternions  as: 


T  1  -  2 (q2  +  q2) 

I 

R  =  |  2 (qxq2  -  q3q4 ) 
[2  {qxq3  +  q2q4 ) 


2(4^2  +  ^3^4) 
1  -  2(ql  +  q 3  ) 
2(?2<?3  -  ?,?„) 


2(‘?1‘?3  -  q2q4) j 
2 (<?293  -  9,?4)  | 

l-2(9l2  +  92)J 


and  vice-versa: 
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^2  =■ 


4  3  =' 


(^23-^32) 

4^4 

(^31-^3) 

4^4 

(Rl2-R2l) 

4^4 


<?4  = 


Similar  equations  exist  in  case  =  0  ,  which  make  use  of  any  of  the  other  three  quaternions  at  the  denominator. 

Given  the  initial  attitude  of  the  spacecraft,  expressed  as  a  quaternion,  the  attitude  at  subsequent  times  can  be 
tracked,  knowing  the  angular  velocity  in  body  axes  from  the  Euler  equations,  by  integrating  the  following  differential 
equations  of  motion  for  the  quaternions: 


with:  n  =1 


d  q 

1 

— 

=  n  q 

dt 

2 

r  0 

I 

*>3 

-®2 

1 

|  -^3 

0 

0 

1 

^2  | 

-2 

"®3 

0  1 

I 

l_-®i 

~-®2 

"®3 

1 

0  J 

2.4.  Validation 

The  diagram  in  Fig.  11  shows  how  the  integration  of  orbital  and  attitude  dynamics  is  preformed,  with  the 
interdependency  of  each  one  of  the  torques  and  forces.  The  propagation  is  initiated  with  the  six  initial  Keplerian 
elements  to  define  the  initial  position  and  velocity,  and  Euler  angles  (roll,  pitch  and  yaw)  and  angular  rate  for  the 
attitude.  The  integrator  is  the  Runge-Kutta  solver  ODE45  provided  in  MATLAB.  The  following  sections  show  some 
of  the  tests  that  were  done  for  validating  the  algorithm  and  its  inner  functions. 
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Fig.  11.  Propagation  of  orbital  and  attitude  dynamics. 


2.4.1.  Solar  radiation  pressure  force 

To  validate  the  solar  radiation  pressure  force,  we  consider  a  thin  flat  plate,  with  same  reflectivity  properties  on 
both  front  and  back  sides.  The  the  position  vector  of  the  object  and  the  Sun  with  respect  to  Earth  are  fixed,  while  the 
incidence  angle  is  varied  over  the  full  circle  (Fig.  12). 


Fig.  12.  Solar  radiation  pressure  force  result  by  varying  incident  angle 

The  front  side  of  the  plate  is  subject  to  SRP  between  0°-  90°  and  270°  -  360°,  while  in  the  remaining  interval 
between  90°  -  270°  the  back  side  is  illuminated.  Due  to  the  cosine  law,  the  maximum  magnitude  of  the  force  is  at  0° 
and  180°  (plate  perpendicular  to  the  sun  direction)  and  null  at  90°  and  270°  (plate  tangent  to  sun  direction). 
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2.4.2.  Shadow  effect 

The  object  is  placed  in  several  positions  with  respect  to  the  Sun  and  the  Earth,  and  it  is  checked  that  the  returned 
value  for  the  SRP  is  either  the  one  in  lighting  conditions  or  zero  when  inside  the  Earth’s  umbra. 

2.4.3.  Orbital  dynamics 

We  consider  a  flat  plate  whose  characteristics  are  in  Table  1.  The  dynamics  of  the  orbit  is  propagated  starting 
from  initial  Keplerian  elements  in  Table  2,  for  10  orbits.  Orbit  of  the  objects  is  simulated  by  equation  (3.9)  and  set 
the  initial  equation  in  Table  2.  Gravitational  perturbation  of  the  Sun,  Earth  and  Moon  are  considered.  The  resulting 
orbit  in  the  physical  space  is  plotted  in  Fig.  13,  and  the  evolution  of  the  5  slowly-changing  Keplerian  elements  over 
time  is  plotted  in  Fig.  14. 

Table  1.  Object  data  _ 


Object  characteristic 

Value 

Mass  (kg) 

10 

Area  (m2) 

1 

Angular  velocity  [x,  y,  z]  (rad/s) 

[0,0,0] 

Specular  reflection  coefficient 

0.6 

Diffuse  reflection  coefficient 

0.26 

Absorption  coefficient 

0.14 

Conduction  (0.25mm) 

36.59xl06 

Table  2.  Initial  conditions 

Keplerian  element 

Value 

Semi-major  axis  (km) 

42,164 

Mean  anomaly  (degrees) 

270° 

Argument  of  perigee  (degrees) 

90° 

Ascending  node  (degrees) 

60° 

Eccentricity 

0.0001 

Inclination  (degrees) 

40° 
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Debris  Orbit 


Fig.  13.  Trajectory. 
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x  io"4  Eccentricity 


Fig.  14.  Keplerian  elements  over  time. 


2.4.4.  Torque-free  attitude  dynamics 

The  attitude  dynamics  is  validated  by  setting  no  external  torques  and  propagating  in  time.  The  moments  of  inertia 
are  arbitrarily  set  as  Ix  =  1.0,  Iy  =  0.8  and  Iz  =  0.5  kg-m2  and  the  initial  angular  velocity  to  co0  =  [5,1,3]  rad/s.  In 
torque-free  motion,  the  angular  momentum  vector  follows  a  periodic  trajectory  as  expected,  plotted  in  Fig.  15.  In 
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addition,  the  angular  velocity  vector  follow  the  trajectory  given  by  the  intersection  of  the  ellipsoid  of  inertia  and  the 
ellipsoid  of  angular  momentum,  which  are  plotted  in  Fig.  16,  matching  to  the  results  in  [32]. 


Angular  momentum  in  body  frame 


Fig.  15.  Trajectory  of  angular  momentum  in  body  frame. 


Angular  momentum  in  body  frame 


Hy(N.km) 


Hx(N.km) 


Fig.  16.  Ellipsoid  of  energy  and  ellipsoid  of  angular  momentum. 


6 


8 


2.4.5.  Coupled  orbit  and  attitude  dynamics 

The  propagation  of  coupled  attitude  and  orbital  dynamics  is  presented  here,  for  10  days  (Fig.  17,  Fig.  18  and  Fig. 
19),  100  days  (Fig.  20,  Fig.  21  and  Fig.  22),  and  1  year  or  365  days  (Fig.  23,  Fig.  24  and  Fig.  25).  Inclination, 
eccentricity  and  longitude  of  ascending  node  for  365  days  change  only  slightly.  The  longitude  of  ascending  node 
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decreases  as  inclination  increases  up  about  0.5°.  The  eccentricity  initially  increases  to  up  to  0.001  and  then  it  goes 
back  to  zero  after  the  full  year.  This  behavior  is  expected  and  it  is  due  to  the  (apparent)  rotation  of  the  sun  around  the 
Earth.  365  days.  The  plots  of  roll,  pitch  and  yaw  Euler  angles  show  tumbling  on  all  axes. 


Debris  Orbit 


Fig.  17. 10  days,  orbit. 
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Eccentricity 


longitude  of  ascending  node 


argument  of  perigee 
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semi-major  axis 


23456789 
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Fig.  18. 10  days,  Keplerian  elements. 
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Roll  (Nadir) 


Pitch  (Crosstrack) 
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Yaw  (RAM) 


Fig.  19. 10  days,  Euler  angles. 
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Debris  Orbit 


Fig.  20. 100  days,  orbit. 
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Fig.  21. 100  days,  Keplerian  elements. 
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Yaw  (RAM) 


Fig.  22. 100  days,  Euler  angles 


Debris  Orbit 


Fig.  23.  365  days,  orbit. 
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Eccentricity 


longitude  of  ascending  node 


argument  of  perigee 


Time  (s) 


semi-major  axis 


Fig.  24.  365  days,  Keplerian  elements. 
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Fig.  25.  365  days,  Euler  angles. 
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3.  Orbital  and  attitude  dynamics  of  flexible  debris  - 

Bernoulli-Euler  model 


The  results  show  that  the  different  initial  conditions  of  attitude  motion  results  in  significant  changes  in  term  of 
orbit  element  and  the  shape  of  the  MLI  debris  influences  unique  attitude  motion  and  the  volume  entity  of  the 
disturbances.  Thus,  it  is  believed  that  the  deformation  is  one  of  the  main  contributing  factors  to  increase  the  accuracy 
of  orbit  determination  of  flexible  debris. 

In  this  paper,  we  will  use  a  Finite  Element  Method  (FEM)  to  model  the  deformation  of  flexible  debris 
considering  only  in  one  plane  for  this  preliminary  study  and  propagate  the  coupling  of  orbit  and  attitude  dynamics 
over  orbits  in  the  GEO  region  in  3D.  We  will  consider  two  configurations:  flat  sheet  and  folded  flat  plate  and  try  to 
understand  their  evolution  subject  to  environmental  perturbations  due  to  Earth  gravitation,  third  body  gravitation 
from  the  Sun  and  Moon  and  the  solar  radiation  pressure. 

The  work  described  in  this  chapter  was  presented  at  the  64th  International  Astronautical  Congress,  Beijing,  China 
[33]. 


3.1.  Finite  element  analysis 

The  Bernoulli-Euler  theory  [34]  is  used  here  to  investigate  the  deformation  of  debris.  Fig.  26  represents  a  beam 
with  two  nodes.  E  is  Young’s  Modulus,  A  is  cross-section  and  I  is  the  second  moment  of  inertia  and  each  node  has 
three  global  degrees  of  freedom,  two  displacement  in  x,  y  axis  and  one  rotation. 


JC 

Fig.  26.  Bernoulli  beam  element  with  2  nodes 


The  vector  of  displacement  is  given  by: 

U  =[U  1  U2  u3  uA  u5  U6]T 
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where  u  l9  u  2  and  u  3  are  x  and  y  displacement  and  rotation  of  the  1st  node  respectively,  u  4 ,  u  5  and  u  6  are  x  and  y 

displacement  and  rotation  of  the  2nd  node  respectively.  In  this  study,  the  debris  will  have  nine  degrees  of  freedom 
space.  To  simplify  only  the  linear  dynamics  will  be  considered: 

mu +  cu + ku  =  F 


where  u  ,  u  and  u  are  the  vectors  of  generalized  displacement,  velocity  and  acceleration  respectively  : 


r*M 

i  i 

>  ] 

rc/,1 

\U  2 

\ui 

\U  2 

1  u  1 

U  =1  -3  | 

1^4  1 

1^4  1 

\Us 
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L^.J 
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Also,  F  is  the  force  matrix  which  consists  of  forces  along  the  x  and  y  axis  and  M : 

r^i 

F 

I  yl  I 
I  M  ,  I 
F  =  |  '  j 

I  I 
I  /r  I 

|  y2  | 


M  is  the  mass  matrix  standard  of  Bernoulli-Euler  Beam  element  in  local  coordinates: 
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3.2.  Coordinate  transformation 


We  can  transform  from  a  body  frame  of  reference  to  an  inertial  frame 


of  reference  by  using  a  transformation 


matrix: 


r  c  50 

-  5  C  0 

0  0  1 

0  0  0 

0  0  0 

L  0  0  0 


0  0  0] 

0  0  0 

0  0  0 

c  s  0 

-5  C  0 

0  o  1 J 


where 
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„  *2  - e  n„A  .  „  y2~yx 

c  =  cos  6  =  -  and  s  =  sin  0  =  - 

L  L 

Thus,  the  stiffness  matrix  and  mass  matrix  in  inertial  coordinates  are  defined  as: 

k  =  tt  kt  and  m  =  tt  M  T 

Finally,  the  most  common  approach  is  to  define  the  damping  matrix  C,  through  Rayleigh  damping  [35],  which 
assumes  a  proportionality  to  the  mass  m  and  stiffness  k  .  This  equation  is  given  by: 

C  =  ccM  +  J3  K 


where  «  and  p  are  proportional  damping  coefficients 


3.3.  Simulation  analysis 

The  aim  of  this  analysis  is  to  compare  the  orbit  propagation  of  a  debris  considered  as  rigid  body  and  the 
propagation  of  a  debris  including  its  deformation.  According  to  most  observations,  MLI  debris  is  found  in  the  GEO 
region.  We  will  assume  the  initial  orbital  parameters  are  defined  as  in  Table  3. 


Table  3  The  initial  six  Keplerian  elements 


Keplerian  elements 

a  (km) 

41,254 

e 

0.0001 

i  (deg) 

30 

£2  (deg) 

45 

co  (deg) 

14 

M  (deg) 

210 

The  initial  attitude  is  defined  with  all  Euler  angles  and  rates  set  to  zero. 

Reflection  and  material  properties  of  multi-layer  insulations  are  based  on  [36].  The  basic  structure  of  MLI  is 
composed  of  a  single  sheet  of  material,  made  of  Polyethylene  terephthalate  (PET),  thickness  of  6  pm,  and  the 
aluminum  coating  is  1000  A  thick  on  both  sides.  This  investigation  assumes  to  be  one  layer  of  PET  to  simplify  the 
model. 

Two  different  forms  of  MLI  are  considered  for  the  simulations:  a  flat  sheet  and  a  flat  sheet  which  is  folded  along 
the  middle  at  an  angle  of  90  degrees  as  in  Fig.  2(a),  dashed  line.  The  dimensions  of  both  configurations  are:  1  meter 
of  width  and  length,  thickness  of  6  pm.  [37,  38].  The  properties  of  the  material  are  shown  in  Table  4  and  this  yields 
an  area  per  mass  ratio  of:  1 19.90  m2/kg. 


Table  4.  Properties  of  PET  material. 


Type 

Mass 

Young’s 

Poisson 

Cs ,  Cd ,  Ca 

material 

Density 

[kg/m3] 

Modulus 

[N/m2] 

’s  ratio 

PET 

1,390 

3.10xl09 

0.38 

0.6  0.26  0.14 
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3.4.  Deformation  result  validation 


To  validate  the  FEM  code  for  deformation  of  flexible  material,  we  define  that  each  feature  consisting  2  elements, 
which  have  3  nodes,  red  points,  and  apply  force,  8.01  le-7  N,  at  the  1st  node,  blue  direction,  following  Fig.  27(a)  and 
we  define  the  boundary  condition  of  degrees  of  freedom  ,nine  degrees.  We  integrate  the  differential  equations  of 
motion  with  ODE45  function  in  MATLAB  to  study  the  deformation  for  3  seconds.  The  results  are  shown  in  Fig.  27 
and  Fig.  28. 

In  Fig.  27(a),  the  1st  node  displacement  moves  to  the  right  side  while  the  3rd  node  move  downwards.  The  1st  node 
of  the  flat  plate  in  Fig.  28(a)  moves  down  and  results  pulls  with  it  the  2nd  and  3nd  nodes  due  to  tension  forces.  These 
results  are  used  to  validate  the  simple  model  of  the  natural  displacements  of  low  weight  and  thin  material.  Due  to  the 
deformation  in  Fig.  27(b)  and  Fig.  28(b),  the  center  of  mass  moves  from  its  initial  position.  As  a  result,  it  affects  the 
attitude  motion  of  the  debris  as  shown  in  the  next  section. 


Deformation 


Center  of  mass  change 


0  12 


0  115 


E, 

»  0.11  - 


o  Os 
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b) 


x  axis  (m) 


Fig.  27  a)  Deformation  of  folded  flat  plate  at  0,  1,  2  and  3  s  by  force  8.011e-7  N  in  the  first  node  in  right 
direction  b)  The  change  of  center  of  mass  of  folded  flat  plate  at  0, 1, 2  and  3  s 
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Fig.  28  a)  Deformation  of  flat  plate  in  0.2,  0.4  and  0.6  s  due  to  a  force  of  8.011e-7N  in  the  first  node  in  -y 
direction  b)  The  change  of  center  of  mass  of  flat  plate 
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3.5.  Results 

The  Earth  gravitation,  third  body  from  the  Sun  and  Moon  and  SRP  radiation  pressure  have  been  applied  to  four 
different  configurations  characteristics,  rigid  flat  sheet,  flexible  flat  sheet,  rigid  folded  plate  and  flexible  folded  plate. 
The  evolution  of  eccentricity  and  inclination  over  3  days  are  shown  in  Fig.  29  and  Fig.  30  while  a  long  term 
propagation  for  over  150  days  are  shown  in  Fig.  3 land  Fig.  32.  In  Fig.  29  (a),  eccentricity  of  rigid  body  grows  up 
consistently  while  those  of  flexible  body  are  different  that  fluctuate  in  same  peak  of  range  and  inclination  and  it  is 
very  obvious  in  Fig.  29(b)  to  see  inconsistency  of  flexible  debris  of  both  initial  geometry  of  flat  plate  and  folded 
plate.  In  Fig.  30(a),  they  have  similar  trend  to  increase  but  the  flexible  debris  are  below  that  of  the  rigid  body  in  Fig. 
30(b).  To  monitor  in  the  long  period,  the  mean  eccentricity  of  the  orbit  of  the  rigid  body  increases  (Fig.  31)  while  in 
the  case  of  flexible  debris  do  not  rise  significantly.  The  mean  inclination  for  the  rigid  debris  in  Fig.  32  increases 
while  for  flexible  debris  it  oscillates  around  a  mean  value.  The  reason  of  this  can  be  found  in  the  change  of  the  cross- 
sectional  area  due  to  the  deformation  of  the  sheet,  effectively  changing  exposed  area  to  solar  radiation  pressure  and 
sun  angles. 


Fig.  29.  a)  Eccentricity  of  flat  plate  and  folded  up  flat  plate  in  both  rigid  and  flexible  cases  b)  Magnify 
eccentricity  period  of  1.42  -  1.44  days 


Fig.  30.  a)  Inclination  of  flat  plate  and  folded  up  flat  plate  in  both  rigid  and  flexible  cases  b)  Magnify 
inclination  period  of  1.42-1.43  days 
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Fig.  31.  a)  Eccentricity  of  flat  plate  in  both  rigid  and  flexible  cases  in  150  days  b)  Eccentricity  of  folded  flat 
plate  in  both  rigid  and  flexible  cases  in  150  days 


Time[days) 


Time(days) 


Fig.  32.  a)  Inclination  of  flat  plate  in  both  rigid  and  flexible  cases  in  150  days  b)  Inclination  of  folded  flat  plate 
in  both  rigid  and  flexible  debris  in  150  days 


In  Fig.  33  and  Fig.  34,  the  attitude  of  the  flexible  debris  in  both  configurations  appears  to  rotate  significantly 
slower  than  for  the  rigid  debris  for  which  the  flat  configuration  presents  the  fastest  displacements.  As  a  result  of 
deformation,  the  center  of  mass  and  moments  of  inertia  are  always  changing,  leading  to  continuously  varying  effects 
from  SRP  and  gravitational  torques.  After  propagating  over  150  days  in  Fig.  35,  Fig.  36,  Fig.  37  and  Fig.  38,  the 
flexible  debris,  shows  periods  of  small  changes  alternating  with  periods  of  high  fluctuations.  This  effect  does  not 
occur  for  the  rigid  debris  which  exhibits  a  more  uniform  behavior  in  time. 
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Fig.  33.  Euler  angles  of  flat  plate  in  case  of  rigid  body  and  flexible  body  in  3  days 
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Fig.  34.  Euler  angles  of  folded  flat  sheet  in  case  of  rigid  body  and  flexible  body  in  3  days 
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Fig.  36.  Euler  angles  of  flexible  flat  plate  in  150  days 
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Fig.  37.  Euler  angles  of  rigid  folded  flat  plate  in  150  days 
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Time(days) 

Fig.  38.  Euler  angles  of  flexible  folded  flat  plate  in  150  days 


3.6.  Conclusion 

This  section  investigated  the  attitude  and  orbital  debris  of  flexible  debris  in  geostationary  region  in  long  term 
period  150  days.  The  flat  plate  and  folded  flat  plate  have  been  simulated  in  rigid  and  flexible  cases  by  using  same 
properties  of  PET. 

In  a  propagation  of  150  days,  some  differences  can  be  noticed  between  the  rigid  and  the  flexible  cases.  The 
deformation  implies  change  in  the  cross-sectional  area  exposed  to  the  sun,  as  well  as  the  center  of  mass  and  pressure, 
and  moments  of  inertia. 

From  this  investigation,  when  an  accurate  model  of  the  flexible  debris  is  provided,  we  can  improve  to  the 
precision  of  prediction  for  the  orbit  determination,  the  orbital  life  time  and  also  investigate  the  material  types  of  the 
debris  from  rotation  rates. 

Finally,  we  underline  that  investigation  did  not  consider  large  deformations,  for  which  non-linear  methods  would 
be  required  to  increase  the  accuracy  of  the  models  for  orbital  and  attitude  prediction.  In  addition,  a  higher-order 
propagator  shall  be  used  to  decrease  the  simulation  time  and  increase  the  accuracy  in  longer  periods. 
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4.  Multi-body  deformation  model 


The  previous  approach  implied  the  strong  assumption  that  deformations  are  small.  While  in  principle  the 
approach  can  be  applied  to  arbitrarily  large  deformations,  the  results  are  less  and  less  accurate,  the  larger  the 
deformations  are.  Because  of  the  thin  membranes  that  are  modelled  in  this  work,  it  is  expected  that  large 
displacements  might  well  be  possible. 

For  this  reason,  we  introduce  here  an  approach  that  aims  at  modelling  a  flexible  thin  sheet  with  full,  non-linear 
equations  of  motion.  In  particular,  the  debris  object  is  modelled  with  a  number  of  lump  masses,  each  representing  a 
section  of  the  thin  membrane,  inter-linked  ideal  tethers.  Suitable  constraints  will  be  set  in  place,  to  allow  the 
deformation  of  the  structure,  as  well  as  modelling  the  internal  forces  that  arise  due  to  the  deformation  (elastic  and 
viscous).  In  addition,  we  will  consider  self- shadowing  effects  on  the  membrane. 

4.1.  Model  and  dynamics 

A  novel  model  for  a  thin,  highly  flexible  MLI-type  membrane  is  proposed  here,  in  terms  of  multi-body  dynamics, 
and  it  is  solved  using  fundamental  Newtonian  mechanics.  Assuming  that  the  membrane  folds  along  one  line,  the 
deformation  can  be  studied  as  a  two-dimensional  problem  (see  Fig.  39  (a)).  The  membrane  is  modelled  as  a  series  of 
lump  masses,  interconnected  with  rigid  rods.  The  masses  act  as  rotational  joints  for  the  plates  and  include  rotational 
spring  and  damper  to  simulate  the  bending  stiffness  of  the  membrane.  The  membrane  is  free  to  move  and  tumble  in 
the  three  dimensional  space,  however  the  folding  lines  are  always  the  same  on  the  membrane,  and  parallel  to  each 
other,  allowing  to  describe  the  deformation  in  a  two-dimensional  plane  perpendicular  to  the  folding  lines.  In  this 
work,  only  three  lump  masses  and  one  folding  line  are  considered.  This  simplified  model,  in  undeformed  (flat)  and 
folded  states  is  illustrated  in  Fig.  39  (c)  and  (c)  respectively.  Representation  in  3D  space  and  reference  angles  are  in 
Fig.  40  and  Fig.  41. 
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y 


b) 


T orsional  spring 


c) 


Fig.  39.  a)  The  membrane  for  2D  deformation  is  modeled  with  series  of  masses  including  torsional  springs  and 
dampers,  b)  Membrane  completely  flat;  c)  Membrane  folded  along  one  folding  line.  The  torsional  spring  in 
the  model  simulates  the  membrane  bending  stiffness. 
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a) 

z  y 


Fig.  40.  a)  the  flexible  model  in  three  dimensions  and  inertial  frame  b)  the  model  in  x-z  view  c)  the  model  in  x- 
y  view. 


Fig.  41.  Angles  defining  the  direction  of  the  imaginary  rods  with  respect  to  an  inertial  frame. 
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In  order  to  write  the  differential  equations  of  motion  for  this  system  of  masses,  we  imagine  to  cut  the  rods  and  we 
introduce  at  each  end  the  axial  reaction  forces  produced  by  the  rods  on  the  masses,  as  in  Fig.  42. 


c) 

Fig.  42.  Free  body  diagram  of  each  lump  mass  a)  1st  free  body  diagram  b)  2nd  free  body  diagram  c)  3rd  free 
body  diagram 


Now,  the  equation  of  motion  for  each  of  the  three  masses  i=  1,  2,  3  is: 

F total  i  +  Ti  =  miXi 

“  FlolaU  +  T,  (2) 

X.  =  - 

m . 


Ftotal  i  is  total  force  acting  on  lump  masses  (  f roral  =  Fext  +  fs.  +  f d  ).  Fext  is  the  external  force  from  conservative 
and  non-conservative  perturbations  (J2,  SRP  and  third  body  from  the  Sun  and  the  Moon),  t.  is  the  tension  force 
generated  by  the  rod,  fsi ds  the  rotational  spring  force  of  ith  mass,  Fd  is  the  rotational  damper  force,  x.  is  total  three 
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axis  acceleration  of  ith  lump  mass  (J2,  SRP  and  third  body  acceleration).  The  rotational  torques  of  the  spring  and  the 
damper  are  found  from  their  corresponding  torques,  acting  at  mass  2,  between  the  first  and  the  second  rod  (Fig.  43). 


b) 

Fig.  43.  Free  body  diagram  of  each  lump  mass  a)  Forces  from  torsional  spring  b)  Forces  from  torsional 
damper 

The  rotational  spring  force  is: 


f  =  k  e 
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EI 

L 


Where  k  is  rotational  spring  constant,  0  is  angle  of  deformation,  e  is  young  modulus,  /  is  the  moment  of 
inertia  of  thin  plate,  l  is  the  length  of  each  rod. 

The  rotational  damping  force  is: 

Fd  =  c  e 

c  =  DF  ^Mks 


Where  c  is  the  rotational  damping  constant,  6  is  angular  velocity  of  the  deformation,  d  f  is  dissipation  factor 
of  material  [39],  m  is  total  mass  of  membrane. 

To  be  able  to  integrate  the  equations  of  motion  (2),  the  rod  tensions  T  is  still  to  be  determined.  To  do  that,  we 
exploit  the  constraint  that  each  rod  is  infinitely  rigid  in  the  axial  direction,  and  thus  (length  of  connecting  rods,  /. ): 

,  .  2  ,  .2,  .2,2 

U,-M  -  x.)  +(yi+i-yi)  +(zi+l-z{)  =  l 


Taking  the  second  differential  of  these  two  equations  we  get: 

2(x.  -  xi+1)(x.  -  xi+l)  +  2 ( y .  -  yM)(yt  -  yi+1)+  2 ( z.  -  zM)(z.  -  zm)  +  2{x.  -  xi+l )2  +  2(y.  -  yi+l )2  +  2(zt  -  zi+l f  =  0 

In  these  equations,  positions  are  velocities  are  known  from  the  state  vector  of  the  system.  The  accelerations  are 
unknown,  however  an  expression  of  them  is  in  Eqs.  (2).  Substituting,  both  tensions  T  appear  in  the  equations,  which 
can  then  be  written  as  a  linear  system: 

rc1  c2irrj  r Aj i 

k  kkkkJ 


Where: 

it  li  li 

C1  =  c1c3(x1  -  x2)( - +  - )  +  c1s3(y1  -  y2)( - + - )  +  -  z2)( - + - ) 

m1  m  2  ml  m  0  ml  m  2 

C2C4  C2S4  S2 

c2  =  --L-t-{xl-x2)--L-*-(yl  -  y2)-—{zx  -  z2) 
m 2  m  2  m2 

C1C3  C1S3  S1 

C3  =  ~—(x2  -  v,)-— (y 2  -  ^3)+  —  (z2  -  z3) 
m  0  m  2  m  3 

11  11  11 

C4  =  c2c4(x 2  -  x3)( - +  - )  +  c2s2(y2  -  y3)( - +  - )  +  s3(z2  -  z3)( - +  - ) 

m  2  m3  m2  m3  m2  m3 

.  ...  F  F  F  !  F  2  p  F 

\  =  ~(xl  -  x2)2  -  (y1  -  y2)2  -  (zt  -  z2)2  -  — t  -  x  2)  +  — ^ —(x  x  -  x  2)  -  ~^{y  x  -  y  2)  +  (y  t  -  y  2)  -  — ^(z  x  -  z  2)  +  x  -  z2)  -  Ft  -  Fd 

m  j  m !  m  j  m  2  mx  m  2 

.  .  .  2  .  .  2  F  F  Fy2  Fy3  f  F 

a2  =  ~(x2  ~  x3)  -  (y2  -  y3)  -  (z2  -  z3)  -  ~^(x  2  -  x  3)  +  ~^~(x  2  -  x  3)  -  — ~(y  2  -  y  3 )  +  —  (y  2  -  y  3)  -  ~^~(z2  -  z3)  +  ~^(z2  -  z  3)  -  F  s  -  F  d 


sx  =  sin^)  s3  =  sin(^3) 

C1  =  cos(^)  C3  =  cos(^3) 

In  which:  and  ,  7?  12  =  -R21  and  R23  =  -R32 

S2  =  sin  (^2 )  =  sin(04) 

c2  =  cos(6>2)  c  4  =  cos((94) 

The  tensions  can  then  be  calculated  by  inverting  the  system: 
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rr,i_rc,  c2r‘rA,i 

kJ~k  c4J  UJ 

Substituting  these  values  in  Eqs.  (2)  allows  to  integrate  them  over  time,  to  determine  the  position  and  velocity  of 
each  lump  mass. 

4.1.1.  Numerical  validation 

To  test  the  spring  and  damper  of  the  model,  it  is  required  to  set  a  non-flat  geometry  as  initial  condition  (triangular 
shape  in  Fig.  44).  The  oscillation  in  Fig.  45  (b)  is  forced  by  the  rotational  spring  only.  The  damper  slowly  reduces  the 
amplitude  of  the  oscillation,  until  it  comes  to  a  complete  stop,  as  in  Fig.  46. 


Fig.  44.  Initial  shape 


Initial  shape 
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motion 

1.2 . . . 


z  axis(m)  x  axis(m)  x  axis(m) 


Fig.  45.  The  oscillation  due  to  deformation  from  non-equilibrium  results  from  torsional  spring  without 
external  force 


motion 


Fig.  46.  Torsional  spring  and  torsional  damper  test  show  that  oscillate  and  back  to  equilibrium  due  to 
torsional  damper. 


4.1.2.  Self  shadowing 

The  planar  shadow  projection,  first  developed  by  Blinn  [40],  allows  shadows  to  be  cast  on  a  planar  surface  as  in 
Fig.  47.  The  point  p  is  the  projection  of  vertex  v  onto  a  plane  P:  n  ■  x  +  d  =  o  due  to  a  location  of  light  source  /  . 
The  plane  P  can  be  described  as 

~  ~  d  +  n  ■  l 

p  =  - - - -(v  -  /) 

n  •  (v  -  /) 
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This  can  be  written  into  a  projection  matrix  (Mv  =  p)  [41]. 


r n.l  +  d  -  l  n 

1  X  X 

-l  n 

x  y 

-l  n 

X  z 

-lxd 

-/  n 

y  x 
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y  y 

-l  n 

y  z 

-l  d 
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-In 
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-l  n 

z  y 

n.l  +  d  -  l  n 

z  z 

-l  d 
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-  n 

y 

-  n 

z 

n.l 

The  self- shadowing  algorithm  has  two  main  steps.  Firstly,  when  two  planes  are  oriented  in  space,  it  is  important 
to  check  which  plane  casts  shadows  on  the  other.  The  vector  from  the  Sun  to  the  centre  of  mass  on  a  panel  will 
determine  which  face  is  closer.  Then,  the  algorithm  determined  if  an  area  of  the  panel  is  exposed  to  SRP  by  checking 
the  intersection  of  area  between  the  shadow  and  the  panel.  This  algorithm  is  calculated  on  each  step  of  integration. 

Fig.  47(a)  presents  a  self-shadowing  simulation.  The  panel  above  (grey)  acts  as  a  shading  panel.  The  bottom 
panel  (brown)  has  a  shadow  cast  upon  its  whole  surface.  The  shaded  area  depends  on  the  location  of  the  light  source 
and  the  model’s  shape.  Fig.  47  (b)  shows  different  shadows  created  by  moving  the  light  source  position  and  Fig.  47 
(c)  by  moving  both  light  position  and  changing  orientation. 


Fig.  47  Self-shadowing  area  simulation  a)  light  source  above  the  model  and  full  shadow  on  the  second  plane  b) 
shadow  area  after  moving  light  source  position  c)  shadow  area  after  rotating. 
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4.2.  Average  solar  radiation  pressure  for  rigid  debris 

The  attitude  of  a  piece  of  debris  will  alter  the  effects  of  solar  radiation  pressure  through  Eq.  (1).  For  sake  of 
comparison,  it  is  beneficial  to  average  the  SRP  force  over  the  possible  tumbling  motion  and  use  this  averaged  force 
to  propagate  over  time  in  the  case  of  the  rigid  debris.  To  reduce  computational  cost  we  decouple  attitude  and  orbital 
dynamics.  The  results  can  then  be  averaged  over  one  orbit  to  determine  the  averaged  dynamics  and  compare  it  with 
the  dynamics  of  the  deformed  model. 

To  model  a  tumbling  piece  of  debris,  it  is  assumed  that  it  can  adopt  any  orientation  in  the  inertial  space. 
Therefore,  the  expected  value  of  the  force  is  equal  to  the  average  force  over  all  possible  orientations.  The  average 
force  is  obtained  by  integrating  over  the  solar  latitude  and  longitude  Eq.  (1)  and  set  s  =  [1,  0,  0]. 

1  n  ln  - 

F  =  - f  f  F  d  A  d  S 

avg  I  J  rad  ,  j  s  s 

0  0 


And  the  equivalent  area  is  expressed  as: 


A 


eq 


F 

avg 


Psp(R) 


Where  psp(R) 


is  solar  pressure  per  unit  area.  Then,  the  average  SRP  can  be  calculated  as: 


F 


AVG 


=  -  A  P  (R ) 

eq  SP  v  ' 


/ 

||x. 


4.3.  Simulation  analysis 

Two  different  kinds  of  MLI  [42]  are  selected  for  analysis:  PET®  and  Kapton®.  PET®  is  coated  aluminium  on  both 
sides  while  Kapton®  is  coated  only  on  one  side.  Their  properties  are  listed  in  Table  5.  The  initial  geometry  is  flat  for 
both  rigid  and  deformed  debris.  All  objects  start  with  the  same  set  of  Keplerian  elements,  which  have  a  semi -major 
axis  41254.0  km  eccentricity  of  0.0001,  inclination  of  5.0  degrees,  argument  of  perigee  of  9.0  degrees,  argument  of 
perigee  of  9.0  degree,  mean  anomaly  270.0  degrees.  Initial  Euler  angles  of  both  planes  are  chosen  from  0,  90  and  - 
180  degrees  (3-1-3  sequence)  and  they  allow  to  orient  the  membrane  at  no  particular  direction.  The  initial  angular 
velocity  set  zero  around  all  directions. 

The  orbit  and  attitude  evolution  of  both  debris  models  has  been  investigated  under  various  perturbations.  J2 
perturbation  is  labelled  with  subscript  “J”,  Sun  and  moon  third  body  gravitational  perturbations  are  labelled  with 
subscript  “g”,  while  the  SRP  force  is  labelled  with  subscript  “s”.  Self-shadowing  is  denoted  by  subscript  “h”.  The 
rigid  body  and  flexible  body  have  subscript  “r”  and  subscript  “f”  respectively. 

Table  5  Properties  of  PET  and  Kapton  [42]. 
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Material  type 

AMR 

[m2/kg] 

Young’s 

Modulus 

[N/m2] 

Cs ,  Cd ,  Ca 

PET 

coated 

111.11 

8.81xl09 

0.60  0.26  0.14 

Kapton 

coated 

uncoate 

d 

26.30 

26.30 

2.50xl09 

0.60  0.26  0.14 

0.00  0.10  0.90 

4.3.1.  Results  and  discussion 

Fig.  48  shows  the  dynamical  evolution  for  the  PET®  over  12  days  under  the  Earth  oblateness  (J2)  and  solar 
radiation  pressure.  Both  eccentricity  and  inclination  evolutions  of  fPETjS;h  lay  between  fPETjS  and  rPETjS.  Fig.  49  is 
analogous,  but  for  Kapton. 


a) 


Fig.  48  Comparison  in  eccentricity  and  inclination  evolution  of  fPETJS ,  fPETjs, h  and  rPETjs  under  J2  and  solar 
radiation  pressure  a)  evolution  over  12  days  b)  evolution  over  0.5  days 
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Fig.  49.  Comparison  in  eccentricity  and  inclination  evolution  of  fKapjs  ,fKapjS,h  and  rKapjs  under  J2  and  solar 
radiation  pressure  a)  evolution  over  12  days  b)  evolution  over  0.5  day 

To  study  the  details  of  physical  dynamics,  we  restrict  the  time  of  the  investigation  during  0.32-3.25  day  (7.12 
minutes)  as  in  Fig.  50.  The  simulation  shows  a  time-lapse  of  the  membrane  deformation  and  attitude  state,  absolute 
acceleration  and  Euler  angles.  It  can  recognize  that  deformed  model  is  tumbling  and  changed  their  geometries  during 
orbiting  in  Fig.  50(a).  In  this  investigation,  SRP  is  the  dominating  force  to  affect  the  dynamics.  The  sun  vector  can  be 
assumed  to  be  constant  due  to  the  investigation  in  a  short  period  of  time.  The  deformation  and  rotation  are  occurred 
when  the  solar  radiation  force  is  unbalanced  on  the  two  flat  surfaces  representing  the  membrane.  Basically,  an 
unbalanced  SRP  force  on  each  flat  surface  depends  on  two  main  factors:  first,  the  difference  in  incident  angle  on 
both  planes  due  to  the  deformation  and  rotation  of  a  flexible  model.  Second,  the  self-shadowing  leading  to  the 
variations  of  exposed  area. 

To  show  the  effect  of  different  sun  incidence  angles  between  two  planes,  without  considering  the  shadowing, 
fPETjS  has  been  simulated  by  neglected  self- shadowing.  Fig.  50(a)  shows  that  an  object  can  be  deformed  due  to  the 
difference  incident  angles.  The  absolute  accelerations  of  both  planes  (Fig.  50(b))  are  never  zero  due  to  SRP  forces. 

Conversely,  Fig.  51(b)  shows  that  the  absolute  acceleration  in  some  periods  gets  to  zero  or  very  small,  when  self¬ 
shadowing  is  considered.  Self-shadowing  can  occur  in  both  planes  during  the  simulation  depending  on  which  one  is 
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closer  to  the  Sun  as  calculated  by  the  self-shadow  algorithm,  and  the  shadow  of  a  panel  can  cover  some  area  of  the 
one,  or  cover  it  completely.  As  the  results  of  case  without  self-shadowing,  fPETjS  is  therefore  exposed  by  higher  the 
solar  radiation  pressure  than  fPETjS,h •  This  leads  to  the  highest  secular  eccentricity  gain  and  largest  inclination 
variations  than  others  in  Fig.  48. 


Day 


c) 

Fig.  50  Simulation  results  of  PET  under  J2  and  solar  radiation  pressure  without  self-shadowing  in  0.320-0.325 
days  (7.12  minutes)  a)  time-lapse  of  deformation  b)  absolute  acceleration  of  both  planes  c)  Euler  angle 
evolution 
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Fig.  51  Simulation  results  of  PET®  under  J2  and  solar  radiation  pressure  with  self -shadowing  in  0.320-0.325 
days  (7.12  minutes)  a)  time-lapse  of  deformation  b)  absolute  acceleration  of  both  planes  c)  Euler  angle 
evolution 


Comparing  the  investigation  of  Kapton®  under  the  same  circumstances  during  0.32-3.25  day  (7.12  minutes),  its 
reflection  properties  are  non-uniform.  The  front  surface  is  coated  with  aluminium  that  is  highly  reflective,  while  on 
the  back  side  the  reflectivity  is  very  low.  When  we  neglect  self-shadowing,  there  are  2  phenomena  explaining  the 
dynamics.  Firstly,  if  both  reflective  sides  of  the  planes  are  exposed  to  SRP,  the  object  will  accelerate  following  the 
direction  of  the  SRP  forces  (plus  the  gravity)  if  the  difference  incident  angles  of  both  planes  is  not  large  enough  to 
deform  it.  If  the  SRP  acts  on  both  back  sides,  the  object  will  essentially  experience  very  little  disturbance  for  orbt 
and  attitude.  Lastly,  when  the  SRP  acts  on  the  front  side  of  one  plane  and  back  side  of  the  another  one,  it  will 
activate  deformation  and  tumbling.  This  leads  the  dynamics  of  fKapjS,  shown  in  Fig.  52. 

Comparing  with  fKapjS)h,  it  is  clear  that  the  debris  has  higher  potential  to  change  a  shape  and  rapidly  rotate  in  Fig. 
53(a)  and  Fig.  53(c)  than  for  fKapjS,  due  to  self-shadowing.  The  inclination  and  eccentricity  evolutions  in  Fig.  49, 
therefore,  illustrates  that  fKapjS,h  shows  larger  inclination  variations  and  secular  eccentricity  variations  than  fKapjS, 
and  rKapjS. 
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Fig.  52.  Simulation  results  of  Kapton®  under  J2  and  solar  radiation  pressure  without  self-shadowing  in  0.320- 
0.325  days  (7.12  minutes)  a)  time-lapse  of  deformation  b)  absolute  acceleration  of  both  planes  c)  Euler  angle 
evolution 
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a)  b) 


Day 


Fig.  53  Simulation  results  of  Kapton®  under  J2  and  solar  radiation  pressure  with  self -shadowing  in  0.320- 
0.325  days  (7.12  minutes)  a)  time-lapse  of  deformation  b)  absolute  acceleration  of  both  planes  c)  Euler  angle 
evolution 


The  investigation  of  inclination  and  eccentricity  under  J2  and  the  third  body  perturbation  in  Fig.  54(b)  illustrates 
that  both  evolutions  of  rigid  and  flexible  models  are  not  significantly  different  in  orbital  dynamics,  as  they  differ  in 
the  order  of  10"5  and  10'8  respectively. 
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x  10'5 


C)  d) 

Fig.  54  Comparison  inclination  and  eccentricity  evolution  of  rPETg,  fPETg,  rKapg  and  fKapg  under  J2  and 
third  body  perturbations  over  12  days  a)  Inclination  evolution  b)  Difference  in  inclination  between  flexible 
model  and  rigid  model,  c)  Eccentricity  evolution  d)  Difference  in  eccentricity  between  flexible  model  and  rigid 

The  next  set  of  figures  show  the  orbital  evolution  under  all  perturbations  for  the  four  object  models  (fPETjgth , 
fKapjg,h,  rPETjg  and  rKapjg )  and  will  consider  self-shadowing  for  the  deformed  model.  Fig.  55  shows  the  trend  in 
eccentricity  and  inclination.  fPETjgth  shows  the  highest  amplitude  of  inclination  oscillations  and  secular  eccentricity 
gain.  Considering  the  same  material  type  in  the  flexible  model,  fPET  has  uniform  reflection  properties  and  larger 
AMR  (because  of  its  lower  density)  than  that  of  Kapton®  (around  4  times)  leading  to  a  greater  inclination  oscillation 
amplitude  and  a  higher  secular  eccentricity  increase.  Moreover,  it  is  evident  from  all  simulations  that  variations  of 
both  inclination  and  eccentricity  caused  by  the  direct  radiation  pressure  are  dominating  over  gravitational 
perturbations. 

Finally,  for  the  self-shadowing  and  flexible  model  of  both  PET®  and  Kapton®,  we  present  for  reference  Fig.  56 
and  Fig.  57,  which  show  their  unique  attitude  dynamics  over  1  day. 
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Fig.  55  Comparison  in  eccentricity  and  inclination  evolution  over  12  days  of  fPETjg,h ,  fKapjg,h ,  rPETjg  and 
rKapjg  under  J2  ,  third  body  perturbations  and  solar  radiation  pressure. 
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c) 

Fig.  56  Euler  angle  evolution  fPETjg.h  over  ldays  under  J2,  third  body  perturbations  and  the  direct  solar 
radiation  pressure  a)  1st  Euler  angle  b)  2nd  Euler  angle  c)  3rd  Euler  angle 
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c) 

Fig.  57  Euler  angle  evolution  fKapjg,h  over  ldays  under  J2,  third  body  perturbations  and  the  direct  solar 
radiation  pressure  a)  1st  Euler  angle  b)  2nd  Euler  angle  c)  3rd  Euler  angle 

4.4.  Conclusion 

A  new  model  for  high  area  to  mass  ratio  (HAMR)  objects  based  on  multi -body  dynamics  that  is  able  to  couple 
deformation,  attitude  and  orbit  dynamics  in  near  GEO  region  has  been  presented.  This  model  uses  a  full,  non-linear 
set  of  equations  for  modelling  the  dynamics  of  a  membrane,  which  is  discretized  with  lump  masses,  a  spring  and  a 
damper. 

The  results  showed  that  perturbations  from  solar  radiation  pressure  play  a  significant  role  in  the  dynamical 
evolution  of  the  debris.  The  J2  and  lunar- solar  perturbations  have  a  significant  impact  on  the  long-term  orbit 
evolution.  The  large  area-to-mass  ratio  is  also  the  reason  of  oscillations  in  inclination  and  secular  increase  in 
eccentricity.  The  flexible  model  shows  a  complicate  tumbling  motion,  which  cannot  be  captured  with  rigid  models. 
Also,  the  self-shadowing  of  the  deformable  model  affects  the  AMR  value  leading  to  the  unique  rotational  and  orbital 
dynamics.  A  significantly  larger  secular  trend  occurs  when  ignoring  self-shadowing  in  case  of  perfect  reflection 
properties.  The  new  model  shows  the  significant  different  variations  of  inclination  and  high  secular  eccentricity 
including  a  rapid  non-uniform  attitude  motion  under  an  influence  of  direct  solar  radiation  pressure  and  gravitation. 
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5.  Propagation  of  coupled  dynamics 

5.1.  Introduction 

In  order  to  study  the  evolution  of  the  orbit  parameters  of  such  objects  over  time,  and  improve  the  understanding 
of  how  does  the  perturbation  affects  the  orbit  parameters,  a  numeric  integration  of  the  orbit  is  commonly  used  [28, 
29]. 

As  discussed  in  the  introduction,  for  long-term  propagation,  up  to  50  years,  the  average  of  the  AMR  value  was 
considered  in  [21].  The  results  were  a  common  trend  of  1  year  period  oscillation  in  eccentricity  in  GEO.  Semi- 
analytical  propagation  was  applied  and  the  value  of  AMR  was  reduced  with  defined  coefficients  in  order  to  take  into 
account  tumbling,  reflection  and  deformation  effects.  Short  periodic  effects  were  not  considered  and  long  periodic 
effects  were  integrated  over  one  revolution.  However,  the  reliability  of  the  propagation  over  time  does  not  seem  to  be 
accurate  not  only  due  to  the  average  of  the  AMR  value,  but  for  possible  resonances  of  the  short  period  terms  such  as 
the  harmonic  of  the  Earth  or  the  Sun -Earth-Moon  configuration  as  it  is  studied  in  [26] . 

An  average  AMR  value  was  also  used  in  [28]  to  propagate  four  objects  from  the  AIUB  catalogue  during  50  days. 
The  ephemerides  of  the  predicted  orbit  were  compared  with  the  observations  of  ZIMLAT,  ESASDT  and  ISON 
network,  and  a  significant  difference  in  distance  values  was  obtained.  This  study  concludes  that  the  values  of  AMR 
are  function  of  time  and  do  not  show  clear  and  obvious  common  trends.  Therefore  further  study  should  be  carried  out 
in  order  to  determine  the  cause  of  the  change  in  AMR  value. 

The  only  non-average  attitude  study  of  MLI-Type  debris  found  in  the  literature  is  carried  out  by  Frlih  in  [10].  The 
results  lead  to  an  overestimation  of  the  AMR  value  if  the  orbital  evolution  is  interpreted  using  the  assumption  of  a 
stable  attitude  rather  than  a  chaotic  tumbling  caused  by  SRP  torque.  The  orbit  propagation  was  carried  out  with  a 
fully  coupled  orbit  and  attitude  dynamics. 

In  all  the  debris  simulations,  the  initial  conditions  for  propagating  the  orbit  are  either  obtained  from  light  curves, 
orbit  break-up  simulations  or  fragmentation  tests  [43].  This  allows  more  realistic  initial  attitude  conditions  and/or 
velocities. 

In  addition  to  the  estimation  of  the  AMR  due  to  deformation,  several  studies  showed  that  different  modelling  of 
eclipse  leads  to  big  errors  when  propagating  the  orbit  [44,  45].  This  effect  seems  to  greatly  affect  MLI-Type  debris  as 
stated  in  [29],  when  perfectly  cylindrical  shadow  cone  and  umbra/penumbra  model  according  to  Baker  were 
compared. 

Clearly,  the  orbit  prediction  of  HAMR  debris  objects,  and  in  particular  MLI-type  debris,  is  a  novel  and 
challenging  field  of  study  where  methods  are  being  proposed  and  validated  with  observation  data,  with  the  aim  of 
increasing  the  understanding  of  this  new  population  of  debris  in  order  to  consider  mitigation  procedures. 

Observations  from  the  ZIMLAT,  ESASDT  and  ISON  network  from  ESA  showed  an  unstable  value  of  AMR  over 
time,  which  were  attributed  to  possible  attitude  motion  or  deformation  (see  Fig.  58,  from  [28]). 
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Fig.  58.  Distribution  of  area  to  mass  ratio  (from  [28]). 


Although  deformation  has  not  been  studied  in  the  literature,  fully  coupled  orbit  and  attitude  motion  of  MLI  debris 
was  investigated  in  Friih  [29].  Fast  spin  motion  was  found  to  be  present  in  case  of  an  offset  between  the  center  of 
pressure  and  the  center  of  mass  of  the  object,  causing  a  high  solar  radiation  torque.  As  it  can  be  seen  in  Fig.  59  (taken 
from  [29]),  the  frequency  of  spin  is  greater  than  the  evolution  of  the  orbit  parameters  in  geosynchronous  orbit. 


Fig.  59.  Euler  angles  function  of  time  (from  [29]). 

This  offset  may  be  caused  by  different  reflectivity  properties  along  the  surface,  for  instance,  due  to  the  tendency 
of  these  objects  to  curl  up  during  peeling  off  -  observed  in  the  degradation  of  the  surface  of  the  Hubble  telescope  -, 
or  due  to  one  end  folded  up  due  to  deformation  -  this  last  one  was  the  one  considered  in  the  simulation  above. 
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Therefore,  a  problem  arise  when  simulating  the  coupled  orbit  and  deformation  of  MLI  debris  with  fast  spin,  as 
the  frequencies  of  the  spinning  result  in  a  severe  increase  of  computational  cost.  Moreover,  the  nature  of  this  type  of 
debris  requires  dealing  with  non-negligible  perturbations,  and  small  integration  steps  for  the  equations  of  motion. 

Fig.  60  shows  the  inter-dependence  of  the  three  dynamics  for  orbit,  attitude  and  deformation.  Each  dynamics  is 
essentially  a  set  of  differential  equations,  which  are  integrated  starting  from  some  given  initial  conditions.  In 
principle,  these  dynamics  are  integrated  separately,  but  their  integration  relies  on  parameters  and  states  of  the  others, 
therefore  cannot  be  simply  decoupled.  Assume  the  SRP  as  the  only  disturbance  to  the  two -body  orbital  dynamics. 
The  SRP  force  can  be  computed  only  when  the  following  are  known:  the  orbital  position  r  (to  give  the  relative 
position  of  the  sun,  eclipses);  the  attitude  of  the  body  with  respect  to  the  sun  q  and  the  state  of  deformation  of  the 
body  s  (defining  the  effective  exposed  area  and  direction  of  the  force).  The  attitude  states  q  and  the  deformation 
states  s  can  be  obtained  integrating  their  respective  dynamics.  However  the  attitude  dynamics  require  knowledge  of 
the  deformation  state  (for  computing  the  inertia  matrix)  and  the  SRP  torque;  the  deformation  dynamics  require 
knowledge  of  the  SRP  forces  and  torques  acting  on  the  body,  as  well  as  the  angular  velocity  «  . 

It  is  also  worth  mentioning  that  in  a  diagram  as  in  Fig.  60,  the  attitude  dynamics  assumes  a  rigid  body,  i.e.  the 
inertia  matrix  of  the  body  changes  slowly  with  respect  to  the  spinning  rate,  and  therefore  any  change  in  the  inertia 
matrix  can  be  neglected  in  the  Euler  equations.  This  assumption  might  not  be  true  for  all  types  of  deformable 
orbiting  bodies. 


Fig.  60.  Dependencies  and  interconnections  among  orbital,  attitude  and  deformation  dynamics. 
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5.2.  Methodology 

The  propagation  of  the  orbit  aims  to  calculate  the  evolution  of  the  orbit  over  time,  starting  from  some  known 
initial  conditions.  The  state  vector  of  variables  to  be  propagated  over  time  can  contain  position  and  velocity 
(Cartesian  elements),  or  the  classical  Keplerian  parameters.  The  first  method  is  known  as  Cowell’s  method,  whilst 
the  second  method  is  referred  as  Variation  of  Parameters  (VOP).  The  VOP  method  usually  presents  singularities 
depending  upon  the  parameters  selected  and  the  orbit  to  be  propagated.  For  this  reason,  other  sets  of  elements  were 
defined  and  used  in  the  literature,  such  as  equinoctial  elements,  which  have  their  singularities  in  points  of  the 
parameter  space  that  are  generally  not  relevant  to  ordinary  orbits). 

When  dealing  with  orbital  perturbations,  it  is  common  to  work  with  the  Cartesian  elements,  since  the 
perturbations,  such  as  the  Solar  Radiation  Pressure  or  drag,  depend  upon  the  position  and/or  velocity  and  cannot  be 
directly  computed  using  other  orbital  parameters,  but  instead  the  orbital  parameters  shall  be  converted  into  Cartesian 
elements  beforehand,  adding  computational  load  on  the  algorithm  at  each  time  step. 

The  propagation  of  the  orbit  corresponds  to  a  mathematic  integration  of  the  selected  differential  equations  of  the 
dynamics  with  respect  to  time,  and  thus  three  techniques  can  be  applied:  analytical,  numerical  and  semi-analytical 
integration  depending  whether  the  integration  can  be  solved  directly  or  not. 

Analytical  methods  show  a  clear  advantage  due  to  the  high  computational  speed.  However,  the  perturbations 
added  to  the  main  gravitational  attraction  must  be  modelled  with  a  mathematical  expression  suitable  for  analytical 
integration.  Generally,  this  mathematical  modelling  is  obtained  by  expressing  the  perturbation  with  series  expansions 
and  truncating  up  to  the  desired  order.  The  truncation  usually  causes  a  fast  degradation  of  the  results  even  with  short 
period  of  propagation.  Most  analytical  techniques  make  use  of  the  variation  of  parameters  proposed  by  Lagrange, 
which  is  only  valid  for  conservative  forces  such  as  the  harmonics  of  the  gravity  field. 

Semi-analytical  methods  are  partly  analytical  and  partly  numerical.  They  exploit  the  advantages  of  both  methods, 
i.e.  fast  computation  of  analytic  methods  and  accurate  expressions  for  perturbation  of  numeric  methods.  The 
propagation  of  the  orbit  is  split  in  two  parts,  an  analytical  integration  for  the  two -body-problem,  and  a  numeric 
integration  for  the  perturbation  dynamics.  At  pre-defined  time  steps,  the  two  dynamics  are  coupled  through  a 
correction. 

This  method  is  efficient  only  if  the  perturbations  are  small  when  compared  to  the  gravitational  attraction,  since 
this  method  relies  on  the  fact  that  the  perturbations  are  evaluated  at  a  lower  frequency,  and  assumed  constant 
throughout  a  time  interval.  If  this  does  not  happen,  then  the  correction  should  be  carried  out  at  a  high  frequency,  and 
the  advantage  of  this  method  is  lost,  and  instead  a  fully  coupled  algorithm  would  be  more  efficient. 

An  example  of  a  semi-analytical  method  is  the  Encke’s  method.  Encke’s  method  splits  the  acceleration 
experienced  by  the  object  in  two  contributions,  one  due  to  the  exact  two -body-problem  dynamics  and  one 
considering  all  other  perturbations: 

r  =  -  —r  +  ap  (r,r ,?,...)  (3) 
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,  the  two-body-problem  can  be  integrated  analytically,  giving  an  accurate  solution  of 

the  unperturbed  (osculating)  orbit  within  a  small  time  interval,  as  seen  in  Fig.  61,  where  p  is  the  radius  of  the 
osculating  orbit,  r  is  the  radius  of  the  actual  orbit,  and  Sr  is  the  difference  in  position  between  the  two  orbits,  i.e. 

r  =  p  +  Sr  . 


initial  position 


osculating  oiN 


Fig.  61.  Definition  of  position  vectors. 


The  osculating  orbit  is  the  exact  solution  of  the  two -body-problem  for  a  given  initial  condition.  It  can  be  shown 
[32],  that  developing  equation  (3)  leads  to: 

JU  JU  ( p3  3  •  , .. 

Sr  =  -  —Sr - -|  - -1  |r  +  ap(r,r,f,...)  (4) 

P  P  \  r  ) 

The  procedure  is  to  firstly  integrate  analytically  the  two-body-problem  to  obtain  the  osculating  orbit  p  (*)  with 

the  initial  state  values.  Then,  Eq.  (4)  is  integrated  to  obtain  £r(f)  until  \sr\  surpasses  a  certain  tolerance.  At  this 

point,  the  algorithm  is  reinitialized  again  setting  the  new  osculating  orbit  as  the  final  real  orbit  obtained  previously. 
The  osculating  orbit  and  the  real  orbit  coincide  at  the  time  of  correction,  when  Sr  =  Sr  =  o  .  The  tolerance  can  be  set 
as: 

krl  |^r| 

H  lrl 

Numerical  methods  are  widely  applied  to  the  propagation  of  the  Cartesian  state  vector:  Cowell’s  method  and 
others.  The  main  advantage  of  this  method  is  that  the  perturbations  are  easily  written  as  function  of  the  state  vector 
variables  -  position  and  velocity  -  and  therefore  they  are  directly  added  to  the  acceleration  term  and  integrated. 

Numerical  integration  methods  can  also  be  applied  to  the  variation  of  parameters,  when  the  perturbations  are 
considered  as  non-conservative  forces,  making  use  of  the  Gauss  equations. 

An  adaptation  of  the  Encke’s  method  has  also  been  proposed  [46]  in  order  to  add  the  known  attitude  dependence 
on  the  orbit  dynamics  (Fig.  62).  The  method  consists  in  integrating  both  the  osculating  and  the  perturbed  orbit 
through  numeric  integration.  However,  the  osculating  orbit  may  not  be  restricted  to  the  two  body  problem,  but  can 
include  orbit  perturbation  terms  -  such  as  higher  harmonics  of  the  gravity  field  -  and  it  is  integrated  with  a  large  time 
step,  whereas  the  perturbed  orbit  takes  into  account  the  other  perturbations  related  to  the  spacecraft  attitude,  and  it  is 
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integrated  with  a  small  time  step.  As  in  Encke’s  method,  a  tolerance  is  introduced  to  reinitialize  the  algorithm  to  a 
new  osculating  orbit. 


Orbit  Integration  Step 


Encke  Correction  Steps 


Orbit 

Update 


Axk+i  = 0 


1  qk 


Time 


9k+1 1 

Updated  Attitude 


Fig.  62.  Integration  scheme. 


5.3.  Numerical  integration 

Numerical  integration  of  the  orbit  is  performed  reducing  the  second-order  differential  equations  of  motion  to  a 
system  of  multiple  first-order  differential  equations  of  the  form  y  =  f  ( y ,  t ) ,  where  y  is  the  state  vector.  Given  the 
initial  conditions  y  0  at  time  t  =  t0 ,  the  problem  is  known  as  initial  value  problem  (I VP). 

Numerical  methods  for  first  order  differential  equations  discretizing  the  time-span  and  determining  the  value  of 
the  state  vector  y  ( t )  at  each  time  instant  by  computing  the  values  of  the  derivative  state  vector  y  at  suitable  time- 

nodes:  these  can  be  previous  points  for  explicit  methods  and  following  points  for  implicit  methods.  Depending  on 
how  many  previous  integration  steps  are  used  within  the  computation  of  the  next  step,  the  methods  are  classified  as 
single-step  methods  or  multi-step  methods  [47]. 

Multi-step  methods  use  information  from  more  than  one  previous  integration  step  in  order  to  predict  and  correct 
the  values  for  the  next  step.  Most  of  the  multi-step  methods  are  based  on  polynomial  interpolation,  and  this  is  the 
main  reason  that  is  restricting  theirs  applications.  Note  that  multi-step  methods  are  not  only  restricted  to  explicit 
methods. 

The  most  common  single-step  methods  are  the  Runge-Kutta  family.  These  methods  are  easily  explained  though 
the  so-called  Butcher  tableau.  Consider  a  Runge-Kutta  method  of  s  number  of  stages,  which  means  that  within  a  time 
step,  there  are  s  number  of  evaluations  of  the  differential  equations  of  motion.  Define  the  node  vector  a  [s  x  1]  as  the 
vector  containing  the  times  at  which  the  dynamics  is  evaluated  within  each  integration  step.  Now  define  the  coupling 
matrix  b  [s  x  s]  as  the  matrix  used  to  compute  the  next  state  vector  within  the  step  as  function  of  the  previous  or  next 
derivative  state  vector.  Finally,  the  weight  vector  c  [1  x  s]  is  the  average  weight  for  each  evaluation.  Using  these 
vectors,  the  Butcher  tableau  is  constructed  as: 

a  b 
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For  example,  for  the  first-order  Euler  method  the  Butcher  tableau  is: 


Higher  order  Runge-Kutta  methods  have  been  formulated,  the  most  common  the  4th  order  “classic  Runge-Kutta 
method”  has  the  following  Butcher  tableau: 


0 

0 

0 

0 

0 

1/2 

1/2 

0 

0 

0 

l/l 

0 

1/2 

0 

0 

1 

0 

0 

1 

0 

1/6 

1/3 

1/3 

1/6 

Note  that  for  implicit  methods  matrix  b  may  have  nonzero  values  on  the  diagonal  or  on  the  upper  triangle.  To 
solve  implicit  methods,  iterations  must  be  carried  out  to  convergence. 

A  really  interesting  family  of  Runge-Kutta  methods  is  the  embedded  Runge-Kutta  family.  These  methods,  also 
called  adaptive  methods,  estimate  the  truncation  error  by  means  of  computing  the  difference  of  the  values  obtained 
from  the  zth  and  (/+l)th  solution.  Once  the  truncation  error  is  estimated,  the  step  size  is  adjusted  according  to  a  semi- 
empirical  formula.  The  same  is  used  to  decide  whether  the  step  shall  be  accepted  or  rejected  and  recomputed  with  a 
new  (smaller)  step  size. 

One  of  the  most  common  adaptive  methods  is  the  Dormand-Prince  4(5),  used  in  MATLAB’s  ode45  solver,  with 
the  following  Butcher  table: 

Table  6.  Dormand-Prince  4(5)  Butcher  table. 


0 

0 

0 

0 

0 

0 

0 

0 

1/5 

1/5 

0 

0 

0 

0 

0 

0 

3/10 

3/40 

9/40 

0 

0 

0 

0 

0 

4/5 

4/45 

-56/15 

32/9 

0 

0 

0 

0 

8/9 

19372/6561 

-25360/2187 

64448/6561 

-212/729 

0 

0 

0 

1 

9017/3168 

-355/33 

46732/5247 

49/176 

-5103/18656 

0 

0 

1 

35/334 

0 

500/1113 

125/192 

-2187/6784 

11/84 

0 

5179/57600 

0 

7571/16695 

393/640 

-92097/339200 

187/2100 

1/40 

35/384 

0 

500/1113 

125/192 

-2187/6784 

11/84 

0 

Note  that  there  are  two  vector  c’s:  namely  one  is  used  in  a  4th  order  step  ( c4 )  and  one  in  a  5th  order  step  ( c5 ).  The 
two  steps  are  performed  simultaneously,  and  their  result  is  compared.  This  comparison  provides  an  estimation  of  the 
truncation  error  eT  performed  by  the  algorithm,  which  is  computed  as  follows: 

=  \ht  (c4  -c5)| 


where  h  is  the  (current)  step  size,  f  is  the  matrix  containing  the  derivatives  of  the  state  vector  evaluated  at  each  time 
inside  the  step,  at  points  specified  by  a.  Finally,  a  correction  for  the  step  size  can  be  computed  as  follows: 
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where  p  is  the  order  of  the  Runge-Kutta  method,  and  eT  allowed  is  the  allowed  error  for  each  step.  In  order  to  avoid  too 

many  rejected  steps,  an  additional  parameter  z  is  added.  Common  values  around  0.85  are  widely  used  [48].  Note  that 
the  previous  expression  is  based  on  the  fact  that  the  local  truncation  error  of  the  Runge-Kutta  methods  is  one  order 
higher  than  the  order  of  the  method. 

Other  algorithms  for  controlling  the  error  of  Runge-Kutta  methods  have  been  proposed,  using  two  steps  of  the 
same  order,  but  of  different  sizes,  for  estimating  the  truncation  error.  However,  these  methods  have  to  compute  a 
higher  number  of  function  evaluations  per  step,  because  the  two  step  sizes  do  not  share  the  evaluation  of  the 
dynamics  at  the  same  time.  Instead,  the  embedded  Runge-Kutta  4(5)  method  discussed  above  requires  5  evaluations 
for  c 5 ,  out  of  which  4  are  shared  with  c4 ,  and  therefore  only  5  evaluations  are  necessary  at  each  step. 

5.4.  Algorithm 

An  adaptive  step,  high-order  Runge-Kutta-Fehlberg  7(8)  integration  method  will  be  implemented  and  validated. 
This  type  of  integration  method  has  been  already  used  in  the  literature  [46]  for  integrating  coupled  orbit-attitude 
dynamics.  An  adaptive  method  such  as  the  one  considered  is  a  good  choice  when  dealing  with  unknown  typical 
frequencies  of  the  dynamics,  and  in  particular  the  attitude  one,  since  the  step  size  is  automatically  readjusted  for  the 
given  admitted  relative  and/or  absolute  errors.  Also,  in  case  of  increasing  or  decreasing  spinning  rate  of  a  debris 
object  over  time,  the  integrator  will  re-adapt  the  step  size  for  optimal  performance. 

The  algorithm  has  been  adapted  from  [48].  The  iterative  integration  procedure  is  shown  in  the  data  flow  in  Fig. 
63.  Essentially  the  differential  equations  of  motion  are  evaluated  at  13  different  times  within  the  given  time  step 
(according  to  the  Runge-Kutta-Fehlberg  7(8)  integration  method),  the  truncation  error  eT  is  estimated  and  compared 

with  the  allowed  error  eTallowed  .  If  the  former  is  smaller  than  the  latter,  then  the  integration  step  is  kept,  and  the  time 
updated,  else  it  is  discarded.  In  both  cases,  the  step  size  is  updated  based  on  ,  and  the  procedure  is  repeated,  until 
the  final  time  is  reached. 
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Fig.  63.  Flow  diagram  for  integration  scheme. 

The  relative  tolerance  RelTol  and  absolute  tolerance  AbsTol  are  the  parameters  used  to  determine  the  accuracy  of 
the  integration,  and  they  are  used  to  compute  the  allowed  truncation  error  eT  allowed  at  each  integration  step,  as 
follows: 


eT  ,allowed 


RelTol  |ym 
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where  ymax  is  the  maximum  value  of  the  state  vector.  However,  to  avoid  numerical  problems  with  relatively  small 

values  of  the  state  vector  -  for  instance  when  determining  the  attitude  in  case  the  Euler  angles  are  nearly  zero  -  a 
minimum  error  will  be  enforced  by  using  the  absolute  tolerance: 

eT  allowed  =  m  1 n  (  ^  ^  ^  ^  |  J  m  ax  |  >  AbsTol^j 

The  step  size  update  is  obtained  with  the  following  equation: 


{eT  +  ePS  ) 

The  small  number  eps  is  added  to  the  truncation  error  since  this  takes  into  account  the  round-off  error  of  the 
calculator.  The  factor  z  is  introduced  in  order  to  tune  the  integrator  and  avoid  unnecessary  further  rejected  steps. 
Note  that  it  must  be  z  <  1  such  that  the  step  size  is  reduced  in  case  of  a  rejected  step.  It  is  set  z  =  0.85  . 

In  addition,  a  minimum  step  size  is  enforced  in  order  to  avoid  step  sizes  becoming  to  small,  which  would  lead  to 
a  sensible  round-off  error.  The  minimum  step  size  is  usually  related  with  the  current  integration  time: 

hmin  =16  eps  (?) 

where  the  function  eps  (?)  is  defined  (in  MATLAB)  as  the  minimum  value  eps  for  which  t  +  eps  =  ?  . 

The  initial  step  size  is  calculated  at  the  beginning  of  the  integration  for  the  first  step,  as  follows: 

h{0)  =  0.5  (RelTol\\y0\\)/S  (5) 

where  the  factor  0.5  again  was  added  to  avoid  o veres timation  of  the  initial  step  size. 

5.5.  Validation 

In  order  to  test  the  implemented  algorithm,  the  integration  of  a  set  of  differential  equations,  for  which  an 
analytical  solution  is  known,  is  considered.  In  this  way,  the  numerical  error  of  the  integrator  can  be  effectively 
measured.  The  error  is  defined  as  the  maximum  difference  between  the  numeric  solution  and  the  analytic  solution  at 
any  time  step  during  the  integration. 

The  proposed  set  of  differential  equations  is  the  following: 

.  ro  n 

"l-  of 

with  initial  condition  y  0  =  [o  i]r  ,  for  which  the  analytical  solution  is  a  unit-amplitude  cosinusoid.  This  is  integrated 
for  100  periods  (final  time  tf  =  200;z-  ).  The  proposed  integration  method  ( rkf78 )  is  compared  in  terms  of  error  and 

function  evaluation  required  against  two  other  numerical  methods:  an  own  implementation  of  a  Runge-Kutta- 
Fehlberg  4(5)  integrator  ( rkf45 )  [49],  and  the  MATLAB  ode45.  These  results  are  presented  in  Table  7  to  Table  10, 
for  different  values  of  AbsTol  and  RelTol.  NFEVAL  is  the  number  of  evaluations  of  the  differential  equations  of 
dynamics,  and  we  assume  here  that  most  of  the  algorithmic  cost  comes  from  evaluating  these  equations  in  a  real 
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case,  and  therefore  the  lower  the  number,  the  better  the  performance  of  the  numerical  integration,  for  the  same 
accuracy.  ERROR  is,  as  mentioned  before,  the  maximum  absolute  difference  between  the  analytical  and  the 
numerical  solution  over  the  integration  time  span.  NSTEP_ACC  is  the  number  of  accepted  integration  steps,  and 
NSTEP_REJ  is  the  number  of  rejected  one  (i.e.  step  was  recomputed  with  a  smaller  step  size).  Finally,  as  a 
reference,  the  computational  time  is  given  (in  seconds).  While  this  value  depends  on  the  machine  used  to  run  the 
algorithms,  it  is  possible  to  compare  the  overall  speed  of  the  three  approaches. 


Table  7.  AbsTol  =  RelTol  =  1.00E-6 

ode45 

rkf45 

rkf78 

NFEVAL 

19357 

20406 

7813 

NSTEP_ACC 

3226 

3187 

601 

NSTEP_REJ 

0 

214 

0 

ERROR 

5.92E-05 

3.25E-04 

4.50E-05 

Time,  s 

1.67 

1.37 

0.452 

Table  8.  AbsTol  =  RelTol  =  1.00E-8 

ode45 

rkf45 

rkf78 

NFEVAL 

49387 

47748 

14313 

NSTEP_ACC 

8231 

7957 

1101 

NSTEP_REJ 

0 

1 

0 

ERROR 

5.60E-07 

3.64E-06 

1.62E-07 

Time,  s 

3.98 

3.46 

0.796 

Table  9.  AbsTol  =  RelTol  =  1.00E-10 


ode45 

rkf45 

rkf78 

NFEVAL 

124705 

120150 

25441 

NSTEP_ACC 

20784 

20024 

1957 

NSTEP_REJ 

0 

1 

0 

ERROR 

5.52E-08 

3.60E-08 

8.59E-10 

Time,  s 

8.97 

10.963 

1.48 

Table  10.  AbsTol  =  RelTol  =  1.00E-12 

ode45 

rkf45 

rkf78 

NFEVAL 

313831 

301896 

45279 

NSTEP_ACC 

52305 

50315 

3483 

NSTEP_REJ 

0 

1 

0 

ERROR 

5.50E-11 

3.59E-10 

4.75E-12 

Time,  s 

24.01 

46.118 

2.621 

It  can  be  seen  from  the  tables  above  that  rkf78  outperforms  the  other  two  4(5)  methods  both  in  terms  of  function 
evaluations  and  error.  The  algorithm  is  considerably  faster  to  run.  ode45  and  rkf45  are  essentially  equivalent  in  terms 
of  number  of  function  evaluations  required,  being  based  on  a  similar  scheme,  even  if  MATLAB’s  own 
implementation  is  optimized  in  terms  of  running  time. 

The  conservative  values  used  for  z  and  the  coefficient  0.5  in  Eq.  (5)  guaranteed  that  no  step  was  rejected  for  any 
of  the  test  cases  in  the  rkf78. 
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6.  Orbit-attitude  propagation  techniques 

6.1.  Introduction 

In  this  chapter,  we  compare  three  different  approaches  for  numerical  propagation  of  coupled  orbital  and  attitude 
dynamics  of  a  body  with  finite  dimensions  (i.e.  no  point  mass).  Essentially  the  idea  is  to  decouple  the  two  dynamics, 
which  might  have  different  characteristic  time  scales,  by  integrating  them  separately.  This  allows  integrating  the 
slower  dynamics  with  a  longer  time  step,  by  achieving  the  same  relative  error  of  the  faster  dynamics,  requiring  a 
smaller  time  step.  The  integration  with  a  longer  time  step  will  allow  saving  computational  time.  Of  course,  the  two 
dynamics  cannot  be  integrated  completely  independently,  due  to  the  fact  that  the  states  of  one  dynamics  are 
parameters  for  the  other  one,  and  vice-versa. 

Note  that  we  are  dealing  with  the  dynamics  of  the  orbit  and  the  attitude,  neglecting  the  dynamics  of  the 
deformation:  we  consider  a  rigid  body.  This  simplifies  the  diagram  in  Fig.  60,  by  removing  the  deformation 
dynamics  and  all  dependencies  from  s  .  This  is  done  for  simplicity.  Whilst  the  methods  and  results  described  in  this 
chapter  are  applied  to  the  coupled  attitude/orbit  problem,  it  is  worth  underlining  that  the  methodologies  proposed  are 
quite  general,  and  can  be  applied  to  any  number  of  coupled  dynamics,  as  long  as  they  have  different  typical  time 
scales. 

In  particular,  the  test  case  adopted  is  a  rigid,  high  area-to-mass  ratio  MLI-type  debris  in  near-geosynchronous 
orbit. 

In  the  following,  a  study  of  three  different  approaches  for  the  coupled  attitude  and  orbit  propagation  problem  will 
be  studied  and  compared  in  terms  of  accuracy  and  computational  efficiency. 

The  first  approach  is  the  integration  of  the  fully  coupled  differential  equations  of  motion  for  orbit  and  attitude. 
Essentially  the  two  dynamics  are  always  evaluated  and  integrated  with  the  same  time  step,  as  if  they  were  one.  In  this 
way,  all  the  states  are  readily  available  at  each  time  step,  without  any  further  calculation.  This  way,  the  time  step  of 
the  integration  will  be  determined  by  the  fastest  dynamics  (namely  the  attitude,  for  a  spinning  body),  and  both 
dynamics  will  be  solved  at  this  time  step.  This  formulation  however,  only  allows  propagation  in  short-term  (few 
days),  due  to  the  computational  cost  required  to  keep  the  integration  error  low.  Full  coupling  with  short-term 
prediction  has  been  studied  by  Frlih  [29] . 

The  second  approach  decouples  the  attitude  and  orbit  dynamics.  The  main  advantage  of  decoupling  the  two 
dynamics  is  the  possibility  to  evaluate  the  slower  dynamics  (the  orbit  in  this  case)  with  a  large  step  size,  therefore 
saving  function  evaluations  and  thus  computational  cost,  as  illustrated  in  Fig.  64.  However,  because  one  of  the  two 
dynamics  is  fast  (namely  the  attitude),  a  small  time-step  size  will  still  be  necessary,  and  therefore  it  is  predicted  that 
even  in  this  case,  the  propagation  would  still  be  limited  to  the  short  term. 
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States 


Fig.  64.  Time  histories  of  the  dynamics  of  orbit  and  attitude. 

The  only  solution  for  long-term  propagation  is  to  eliminate  the  fast  dynamics  of  the  problem.  This  allows  using 
an  overall  larger  step  size,  while  maintaining  the  propagation  accuracy.  Fast  dynamics  can  be  eliminated  with 
averaging  techniques,  when  the  fast  motion  is  periodic.  This  is  certainly  the  case  for  a  spinning  rigid  body. 
Nevertheless,  in  this  case,  while  the  propagation  itself  can  be  very  accurate,  errors  can  still  be  present  due  to  the 
averaging  technique  used.  Averaging  methods  for  long-term  -  up  to  20  years  -  were  used,  for  example,  in  [23]. 

It  is  worth  mentioning  that  currently,  techniques  combining  orbit  prediction  and  observational  data  are  being  used 
for  mid-term  propagation  (~  50  days).  Essentially  objects  are  tracked  with  a  telescope  to  get  an  estimation  of  their 
position  and  velocity,  when  they  are  visible.  These  data  are  then  used  to  propagate  the  orbit  and  predict  the  position 
of  the  next  possible  observation.  This  report  however  is  focused  on  the  propagation  of  the  orbit  only,  assuming  some 
initial  data  is  available  (e.g.  from  a  previous  observation)  and  assuming  that  no  other  observational  data  will  be 
available  during  the  propagation  time  (i.e.  relying  only  on  the  initial  conditions). 

In  the  following  section,  a  prototype  of  the  problem  is  proposed,  and  the  three  different  approaches  (fully- 
coupled,  decoupled  and  averaged)  are  studied  in  terms  of  accuracy. 

6.2.  Formulation  of  the  problem 

For  an  accurate  propagation  in  geosynchronous  orbit,  the  effects  of  third  bodies  (Moon  and  Sun),  the  harmonics 
of  the  Earth,  and  the  perturbation  of  the  SRP  must  be  considered.  For  example,  Friih  [29],  considered  the  harmonics 
of  the  Earth  up  to  sixth  order.  In  the  same  work,  it  is  also  found  that,  regarding  the  attitude  motion,  the  SRP  Torque 
is  by  means  the  most  significant  attitude  perturbation,  while  the  magnetic  torque  seems  to  be  negligible.  The  gravity 
gradient  can  be  significant  depending  on  the  cyclic  differences  of  the  principal  moments  of  inertia  of  the  body,  which 
in  turn  depends  on  the  body  shape. 

Because  the  aim  of  this  chapter  is  to  compare  the  different  propagation  techniques,  as  opposed  to  achieve  a 
precise  modelling  of  the  dynamics  and  the  perturbations,  a  prototype  problem  is  created  in  the  following  way.  First 
of  all,  the  problem  is  planar.  This  allows  reducing  the  number  of  states  without  changing  the  nature  of  the  problem 
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itself.  The  body  is  a  rigid  flat  plate.  The  orbit  dynamics  is  essentially  the  planar  two -body  problem,  perturbed  by  the 
gravitational  attraction  of  the  Sun  and  the  Moon  (here  fixed  in  an  inertial  Earth -centered  frame),  and  the  SRP  force. 
The  attitude  dynamics  is  based  on  the  classic  rotational  equations  of  motion,  and  forced  by  the  SRP  torque. 

The  debris  object  is  modeled  as  a  rigid,  square  flat  plate,  of  i  =  \m  side,  20  microns  thick,  and  the  properties  of 
the  flat  plate  are  as  shown  on  Table  11.  The  area-to-mass  ratio  (AMR)  (inverse  of  areal  density)  is  that  of  PET,  111 
m2/kg.  This  allows  to  compute  the  total  mass  and  the  moment  of  inertia  of  the  object. 

In  order  to  introduce  torques  due  to  SRP,  we  assume  that  the  surface  of  the  flat  plate  has  two  different  parts,  with 
different  reflectivity  properties  (see  Fig.  65).  Most  area  of  the  debris  is  specularly  reflective  (coefficient  of 
reflectivity  1),  while  one  strip  is  fully  absorbing  the  radiation  (coefficient  of  reflectivity  0).  Named  d  the  fraction  of 
total  surface  that  is  non-reflective,  the  equivalent  coefficient  of  reflectivity  of  the  object  is: 

c  r  =  1  -  d 

Different  values  of  d  will  be  used  for  the  The  solar  radiation  pressure  force  can  be  computed  as: 

=  G*f(l  +  O-cos2(0)n 
m 

where  psrp  =wjc  with  solar  constant  wo  =  1 3 68W  /  m 2  at  Earth  distance  and  c  =  3  x  l  o8  m  /  s  the  speed  of  light  in 
vacuum,  e  is  the  angle  between  the  unit  normal  to  the  surface  n  and  the  Sun  direction  (see  Fig.  66).  The  sign  of 
fsrp  must  be  adjusted  to  make  sure  that  it  is  always  pointing  in  the  semi-space  away  from  the  Sun. 


◄ - ► 

Fig.  65.  Flat  plate. 

The  effect  of  the  non-reflective  area  is  also  to  shift  the  center  of  pressure  with  respect  to  the  center  of  mass 
(which  is  in  the  geometric  center),  by: 

r  =  ldl  2 

cmcp  / 

This  shift  generates  a  solar  radiation  pressure  torque  tsrp  =  fsrp  id  /2  . 
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Table  11.  Flat  plate  data. _ 

Object  side,  l  1  m 

Object  thickness  20  microns 

Area,  A  1  nr 

AMR  (area-to-mass  ratio),  A/m  111  m2/kg  (PET) 

Mass  9  g 

Coefficient  of  reflectivity,  reflective  side  1 

Coefficient  of  reflectivity,  non-reflective  side  0 


The  problem  is  then  reduced  to  three  degrees  of  freedom  only:  two  translational  and  one  rotational.  The  states  to 
be  integrated  states  include  these  three  degrees  of  freedom,  and  their  derivatives. 


Fig.  66.  Picture  of  orbital  motion  and  definition  of  attitude  angle  0 . 

Naming  the  ECI  orbital  position  vector  r  ,  and  the  attitude  parameter  e  ,  we  have  that  the  net  acceleration  on  the 
body  is: 

a  =  Mr)  +  a,(r)  +  aM  (r)  +  a^p(r’^) 

This  equation  highlights  the  dependence  of  the  acceleration  due  to  Earth,  Sun  and  Moon  on  the  orbital  position 
r  ,  while  the  SRP  acceleration  is  responsible  for  the  coupling  of  the  dynamics,  in  that  it  depends  on  both  r  and 
attitude  e  .  Similarly,  the  total  torque  acting  on  the  body  is: 

M  =Ms*p(r’0)  +  M*(r’0) 

which  highlights  that  both  terms  -  gravity  gradient  torque  and  SRP  torque  -  depend  on  both  orbital  position  and 
attitude.  Both  attitude  and  orbital  dynamics  are  inter-dependent  upon  each  other,  which  is  representative  of  a  general 
case. 

The  integration  is  performed  in  an  Earth -centered  inertial  (ECI)  reference  frame,  with  Cartesian  coordinates, 
where  the  x-y  plane  is  arbitrarily  the  equatorial  plane. 

The  initial  orbit  is  geosynchronous,  with  initial  state: 
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[42241.08007]  [  0  ] 

rQ  =  j  0  j  km;  vQ  =  j  3.07  1  858028  j  km  /  s 

L  o  j  L  o  j 

=  0 

and  the  equations  are  integrated  for24/z  =  86400s. 

Due  to  the  integration  lasting  1  day  only,  the  position  of  the  Sun  and  the  Moon  are  fixed  in  an  Earth-centered 
inertial  frame.  This  assumption  would  not  be  valid  for  longer  integration  times,  but  it  does  not  affect  the  results 
presented  here  nonetheless.  The  position  of  the  Sun  and  the  Moon  (of  gravitational  constants  1.32712428ell 
kmA3/sA2  and  4902.799  kmA3/sA2)  in  ECI  are  respectively: 

[-149,597,870]  [384,399] 

1  n  I.  I  n  I* 

r5  =  |  0  \km’  rM  =  |  0  \km 

L  o  j  L  0  J 


6.3.  Description  of  the  three  approaches 

6.3.1.  Fully-coupled  approach 

Fully-coupled  dynamics  is  the  conventional  way  to  integrate  different  dynamics  in  case  they  depend  each  other. 
However,  in  case  of  really  fast  spinning  of  MLI-Type  debris,  this  approach  is  limited  to  short  prediction  times  for 
two  reasons:  firstly,  the  computational  cost  is  high,  and  secondly,  the  integrator  error  increases  quickly.  Fig.  67  is  a 
diagram  showing  the  data  flow  of  a  fully-coupled  dynamics  simulation. 


Fig.  67.  Flow  diagram  of  fully-coupled  orbital  and  attitude  dynamics. 
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In  addition,  due  to  the  particular  configuration  of  the  flat  plate  chosen  for  this  simulation,  the  spin  rate  increases 
continuously  due  to  the  SRP  torque,  and  this  causes  the  step  size  to  reduce  along  the  simulation.  Hence,  because  of 
this  increasing  angular  velocity,  the  process  of  averaging  the  attitude  motion  must  be  studied  in  detail  in  order  to 
predict  long-term  behavior  of  this  kind  of  debris. 

6.3.2.  Decoupled  approach 

Decoupled  approaches  for  integrating  dynamics  with  different  timescales  were  found  in  the  literature  of  electric 
circuit  simulation  (e.g.  [50])  for  decoupling  electric  signals  with  different  time  scales.  Essentially,  the  approach  is  to 
use  multi-rate  Runge-Kutta  methods.  These  integrators  are  fixed-step  methods,  and  the  coupling  of  the  signals  is 
made  by  interpolating  or  extrapolating  the  signal  which  is  not  being  propagated  for  a  given  time.  In  a  similar  way,  a 
decoupled  Runge-Kutta  method  is  proposed  and  studied  in  terms  of  accuracy,  and  compared  to  the  fully-coupled 
approach,  which  uses  a  conventional  Runge-Kutta  scheme. 

As  stated  before,  the  core  of  this  formulation  is  in  that  it  allows  evaluating  the  (fast)  attitude  dynamics  with  a 
small  time  step,  and  at  the  same  time  skipping  some  evaluations  of  the  orbital  dynamics  differential  equations,  saving 
computational  cost. 

6.3.2. 1  Formulation 

We  make  use  of  the  Encke’s  method  of  propagation  of  the  osculating  orbit,  and  in  first  instance  consider  that 
over  one  time  step  of  the  orbital  dynamics  ah  ,  the  two-body  problem  is  a  good  approximation,  and  neglect  the 
perturbations  and  the  effect  of  the  change  in  attitude  during  this  time.  This  implies  that  the  perturbations  over  the 
classical  two-body  problem  dynamics  are  very  small  compared  to  the  gravitational  acceleration.  The  time  step  ah 
is  several  times  greater  than  that  of  the  attitude  dynamics,  which  is  integrated  with  step  Ah  .  The  orbit  estimation 
allows  calculating  the  evolution  of  position  and  velocity  after  time  ah  ,  and  at  each  instant  of  time  within  ah 
through  interpolation.  The  idea  is  that,  despite  this  orbital  prediction  is  erroneous,  due  to  the  fact  that  any  variation  is 
attitude  (and  therefore  variation  in  SRP)  is  not  taken  into  account,  on  the  other  hand  these  variations  would  produce  a 
negligible  impact  on  the  single  step.  In  particular,  the  impact  is  negligible  with  respect  to  the  attitude  calculation,  in 
that  the  torques  do  not  change  dramatically  for  minor  variations  in  the  orbital  state. 

This  first  estimation  of  the  orbital  state  is  used  to  propagate  the  attitude  dynamics  along  the  orbital  time  step 
(ah  ),  this  time  using  a  number  of  smaller  steps  Ah  (refer  again  to  Fig.  64).  Despite  the  attitude  propagation  uses  a 
slightly  erroneous  orbital  prediction,  as  discussed  before  it  is  assumed  that  this  error  is  negligible,  and  therefore  the 
attitude  history  obtained  at  this  step  is  final. 

Finally,  the  attitude  time  history  over  ah  is  interpolated  and  used  to  correct  the  position  and  velocity  prediction 
obtained  before. 

The  entire  procedure  is  started  with  the  orbital  and  attitude  initial  conditions,  and  iterates  for  each  step  ah  .  One 
of  these  steps  is  represented  as  a  flow  diagram  in  Fig.  68.  Note  that  as  Ah  <  ah  ,  several  attitude  steps  are  required 
to  cover  one  single  orbital  step.  Therefore,  the  number  of  evaluations  of  the  orbital  (or  slow)  dynamics  is  smaller 
than  that  of  the  attitude  (or  fast)  dynamics.  Finally,  note  that,  because  the  method  is  adaptive,  it  is  possible  that  in 
order  to  not  surpass  the  allowed  truncation  error,  a  h  shall  be  reduced. 
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Single  step  over  ah 


Multiple  steps  Ah 


Single  correction 
step  over  ah 


Fig.  68.  Flow  diagram  of  decoupled  approach. 
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6.3.3.  Averaged  approach 

In  both  previous  approaches  (coupled  and  decoupled  propagation)  the  attitude  has  to  be  propagated  along  time 
with  the  high  frequency.  This  might  lead  to  a  fast  degradation  of  the  accuracy  due  to  the  amount  of  steps.  Averaging 
the  fast  dynamics,  instead,  allows  eliminating  the  fast  frequencies  completely.  While  the  simulation  time  decreases, 
errors  are  introduced  due  to  the  averaging  process,  and  depending  on  the  selected  averaging  method  itself. 

Once  the  averaging  law  is  used,  the  attitude  motion  is  essentially  not  integrated,  and  only  its  net  effect  on  the 
orbit  is  included  in  the  orbital  dynamics,  in  terms  of  accelerations: 

a  =  aE(r)  +  ai(r)+aM  (r)  +  assp  (r,rT) 

Where  it  was  made  explicit  that  the  SRP  term  aSRP  ( r  ,0  ) ,  responsible  for  the  attitude-orbit  coupling,  is  averaged 

over  the  attitude.  Note  that,  in  this  particular  case  of  the  SRP,  not  only  the  area-to-mass  ratio  has  to  be  averaged,  but 
also  the  direction  of  the  solar  radiation  pressure  force. 

Averaging  can  be  obtained  in  different  ways.  In  a  three-dimensional  case  of  randomly  tumbling  object,  and  in 
absence  of  further  data,  it  is  reasonable  to  assume  that  the  object  will  spend  the  same  time  in  each  one  of  all  possible 
attitude  states,  and  therefore  the  resulting  average  exposed  area-to-mass  ratio  is  simply  the  average  of  the  area-to- 
mass  ratio  of  each  possible  attitude  over  time,  and  the  direction  of  the  SRP  force  is  in  average  pointing  in  the 
direction  opposite  to  the  sun.  However,  on  another  assumption,  the  debris  might  be  spinning  around  one  axis,  with 
no  or  small  precession  and  nutation.  In  this  case,  the  direction  of  the  average  SRP  force  might  not  point  away  from 
the  sun,  but  instead  it  has  to  be  evaluated  through  averaging  process. 

For  this  2D  test  case,  the  value  of  the  AMR  will  be  estimated  as  the  average  between  the  maximum  AMR  (flat 
plate  perpendicular  to  sun  line),  and  the  minimum  AMR  (i.e.  zero  when  the  flat  plate  is  aligned  with  the  sun 
direction)  that  the  debris  can  adopt  in  any  angular  orientation.  Therefore,  the  average  SRP  force  is  in  the  sun 
direction  only. 

6.4.  Results 

Three  defferent  test  cases  are  performed,  in  which  the  integration  was  run  using  a  different  value  of  d :  0.01  (1% 
of  object  surface  is  non-reflective),  0.001  (0.1%)  and  0.0005  (0.05%).  Table  12,  Table  13  and  Table  14  show  the 
results  for  each  one  of  these  cases.  Essentially  different  values  of  d  change  the  amount  of  torque  forcing  the  attitude 
motion,  and  therefore  the  typical  frequency  of  the  oscillation  itself. 

Each  of  the  test  case  was  perfomed  in  the  following  way:  first,  the  reference  orbit  and  attitude  time  histories  were 
computed.  This  was  done  by  integrating  the  fully  coupled  equations  of  orbit  and  attitude  with  a  relative  tolerance  of 
10E-16.  The  tables  show  the  number  of  function  evaluations  necessary  to  complete  the  integration,  as  a  reference 
(Full  problem  NFEVAL).  Then,  the  three  different  approaches  were  run  on  the  problem:  fully  coupled,  decoupled 
and  averaged.  For  the  three  approaches,  the  relative  tolerance  RelTol  was  tuned  such  to  reach  a  similar  value  of  the 
final  position  error,  defined  as: 
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error  =  - - 1 

ro 

that  is,  the  difference  between  the  position  at  the  final  time  resulting  from  the  reference  orbit  calculation  and  the 
same  resulting  from  using  each  one  of  the  three  approaches.  The  value  of  ERROR  is  shown  in  the  tables.  In  this  way, 
all  three  approaches  find  the  orbit  with  similar  accuracy,  and  hence  the  number  of  function  evaluations  that  were 
used  can  be  compared.  The  fully  coupled  approach  evaluates  the  orbit  and  attitude  dynamics  simultaneously.  Instead, 
the  decoupled  approach  evaluates  the  orbit  and  the  attitude  dynamics  separately,  and  therefore  the  number  of 
function  evaluations  for  each  one  of  these  (NFEVAL_0  and  NFEVAL_A  respectively).  The  averaged  approach  does 
not  solve  the  attitude  motion,  and  therefore  only  the  orbital  dynamics  is  evaluated. 


Table  12.  d  =  0.01  (1%) 


Full  problem  NFEVAL 

1,192,672 

Coupled 

Decoupled 

Averaged 

RelTol 

5e-14 

1.0E-14 

IE-18 

NFEVAF_A  (c2) 

639,834 

761,605 

/ 

NFEVAF_0  (cl) 

278,837 

92,989 

ERROR 

1.139E-07 

1.250E-07 

1.56E-2 

Table  13.  d  =  0.001  (0.1%) 


Full  problem  NFEVAL 

367,835 

Coupled 

Decoupled 

Averaged 

RelTol 

IE-13 

l.le-14 

IE-18 

NFEVAE_A  (c2) 

183,703 

230,464 

/ 

NFEVAE_0  (cl) 

110,838 

94,315 

ERROR 

1.34e-08 

2.378E-08 

1.56E-2 

Table  14.  d  =  0.0005  (0.05%) _ 

Full  problem  NFEVAL_ 256,412 


Coupled  Decoupled  Averaged 


RelTol 

2E-14 

IE-14 

IE-18 

NFEVAL_A  (c2) 

151,450 

160,095 

/ 

NFEVAL_0  (cl) 

90,142 

94,250 

ERROR 

2.020E-08 

2.371E-08 

1.56E-2 

Focusing  on  the  case  d  =  0.01  first  (Table  12),  the  fully-coupled  approach  requires  639,834  function  evaluations 
to  achieve  a  relative  error  of  IE-7.  The  decoupled  approach  requires  761,605  attitude  equation  evaluations  and 
278,837  orbit  equation  evaluations.  The  fact  that  the  attitude  requires  a  higher  number  of  function  evaluations  is  to 
be  expected,  its  dynamics  is  faster,  and  treated  as  such  by  the  algorithm.  The  decoupled  approach  requires,  overall, 
more  function  evaluations  than  the  fully-coupled  approach,  however  it  is  expected  that  the  fully-coupled  dynamics 
are  more  complex  (and  therefore  more  computationally  demanding)  than  the  attitude  or  orbit  alone.  If,  for  example, 
we  assume  that  the  computational  time  of  a  single  evaluation  of  the  fully  coupled  dynamics  is  unitary,  and  equally 
split  between  calculation  of  the  orbit  and  the  attitude,  then  in  a  first  approximation  we  can  say  that  the  computational 
time  for  each  evaluation  of  each  decouples  dynamics  is  0.5.  In  this  case,  the  computational  effort  of  the  decoupled 
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approach  would  be  520,221  (compared  to  639,834  of  the  fully  coupled  approach),  which  corresponds  to  a  saving  of 
about  18%  of  the  computational  time.  These  figures  are  highly  dependent  on  the  complexity  of  evaluating  the 
attitude  dynamics  with  respect  to  the  orbital  dynamics.  If  evaluating  the  attitude  differential  equations  results  to  be 
much  more  expensive  than  evaluating  the  orbit  equations,  then  the  decoupled  approach  becomes  less  convenient,  or 
even  inconvenient,  with  respect  to  the  fully-coupled  approach. 

Finally,  the  averaged  approach  cannot  reach  the  required  position  error,  stopping  at  1.56E-2  even  if  the  relative 
tolerance  is  set  to  IE- 18.  This  is  certainly  due  to  the  error  introduced  by  not  modelling  the  attitude  dynamics  at  all. 
However,  this  approach  requires  92,989  orbit  function  evaluations,  which  is  only  a  small  fraction  of  what  needed  by 
the  other  two  approaches:  if  high  accuracy  is  not  the  main  requirement,  then  the  averaged  approach  is  without  doubt 
the  most  efficient.  In  addition,  a  different  averaging  could  lead  to  a  smaller  error,  at  no  or  small  additional 
computational  cost. 

The  other  test  cases,  with  different  values  of  d ,  show  a  similar  trend,  however  as  the  typical  frequency  of  the 
attitude  dynamics  lowers,  the  advantage  of  the  decoupled  approach  over  the  fully  coupled  approach  reduces. 
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7.  Conclusions 


This  report  has  presented  two  approaches  to  model  and  integrate  deformable,  highly  reflective  membranes  in 
Earth  orbits.  One  is  based  on  modelling  the  flexibility  of  the  membrane  as  an  Euler-Bernoulli  truss,  the  second 
discretizes  the  membrane  with  lump  masses  connected  by  rods.  The  former  assumes  implicitly  small  displacements, 
which  are  unrealistic  for  large,  thin  membranes  like  MLI  sheets,  while  the  latter  models  arbitrarily  large 
deformations. 

Both  approaches  assume  that  the  membrane  can  fold  along  one  folding  line,  which  remains  fixed  on  the 
membrane.  However,  both  approaches  can  be  extended  by  considering  more  folding  lines,  at  additional 
computational  cost. 

Both  approaches  were  used  to  predict  the  orbit  of  a  typical  GEO  debris,  and  compared  with  other  propagation 
methods  using  rigid-body  approximation,  and  showed  sensible  differences  in  the  evolution  of  the  orbital  parameters, 
particularly  inclination  and  eccentricity. 

The  linear  approach  is  somewhat  intrinsically  limited  due  the  fact  that  large  displacements  cannot  be  correctly 
model  large  membranes,  whose  deformations  are  clearly  not  limited  in  principle.  The  multi -body  approach 
overcomes  this  limitation,  at  the  cos  of  additional  computational  effort. 

Earth  gravity  (up  to  J2),  SRP,  eclipses,  gravitational  attraction  of  sun  and  Moon,  self-shadowing  effects  were  all 
considered  in  the  propagation.  However,  as  it  was  mentioned  in  the  introduction,  because  MLI  membranes  are 
aluminized  and  hence  conductive,  self-induced  currents  are  also  an  important  source  of  internal  forces  that  contribute 
to  the  attitude  dynamics  and  the  deformation  of  the  sheet. 

Despite  all  flexible-body  methods  were  matched  to  equivalent  rigid-body  approximations,  and  the  results  of  the 
orbital  integration  were  compared,  no  validation  of  either  approach  was  done  using  real  space  debris  data  from 
empirical  observations.  Potentially,  light  curves  of  flexible  debris  could  be  estimated  by  considering  the  relative 
position  and  attitude  with  respect  to  the  sun  and  an  observer  on  the  Earth.  Simulations  could  be  carried  out  assuming 
that  a  particular  debris  object  is  a  deformable  MLI  sheet,  and  then  matching  the  simulated  attitude/deformation 
motion  with  the  observed  light  curve.  However,  these  would  depend  on  extremely  precise  initial  conditions,  as  well 
as  characteristics  (size,  thickness,  material)  of  the  debris  sheet  which  are  not  available,  making  this  type  of  validation 
extremely  difficult. 

In  addition,  the  multi-body  approach  makes  use  of  discrete  rotational  springs  and  dampers  to  account  for  the 
bending  stiffness  of  the  membrane.  The  spring  and  damper  constants  were  estimated  through  semi -empirical 
formulas,  based  on  the  data  of  the  materials  of  the  typical  MLI  membranes.  These  estimations,  however,  can  be 
wrong,  because  of  the  extremely  thin  nature  of  the  membrane  itself.  An  experiment  could  be  set  up,  in  order  to  have 
a  better  estimation  of  these  parameters.  This  is  sketched  in  Fig.  69:  an  MLI  sheet  is  held  hanging  on  the  top  wall  of  a 
vacuum  chamber.  The  sheet  is  then  displaced  from  its  equilibrium  position  along  the  vertical,  and  left  swinging  due 
to  gravity.  The  amplitude  and  frequency  of  oscillation  over  time  (monitored  through  a  glass  window)  are  related  to 
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the  elastic  properties  of  the  membrane.  The  vacuum  chamber  ensures  that  the  effect  of  the  air  friction  is  negligible. 
The  spring  and  damper  in  the  multi -body  model  can  then  tuned  such  that  the  response  of  the  model  matches  that  of 
the  real  membrane. 

In  a  second  set  up,  a  beam  of  light  can  light  the  membrane  through  the  glass  window,  and  displacement  of  the 
MLI  membrane  due  to  SRP  can  be  registered,  with  respect  to  the  local  vertical.  According  to  preliminary 
calculations,  this  could  be  in  the  order  of  microns  or  fractions  of  a  micron,  hence  highly  precise  and  sensitive  sensors 
have  to  be  employed,  to  measure  the  displacement  without  contact.  This  variant  of  the  experiment  might  allow  to 
compute  the  effective  SRP  force  on  the  membrane,  once  the  light  flux  through  the  window  is  known. 


Glass  window 


Light  source 


1 


\ 


Fig.  69  Experimental  set  up  for  determination  of  MLI  properties  in  a  vacuum  chamber. 

Finally,  a  numerical  integrator  was  developed  to  integrate  two  coupled  dynamics  with  different  time -scales.  The 
principle  is  that  the  slow  dynamics  can  be  integrated  with  a  longer  time-step  with  respect  to  the  fast  one,  and  its 
states  are  estimated  in  order  to  resolve  the  fast  dynamics  at  each  evaluation.  The  integrator  is  of  general  purpose,  but 
it  was  applied  to  a  flat  plate  with  fast  spinning  dynamics  and  slow  orbital  dynamics. 

It  was  shown  that  a  decoupled  approach  can  have  some  advantages  in  computational  time,  when  integrating 
coupled  attitude  and  orbital  dynamics,  but  with  the  present  algorithm,  the  advantage  is  dependent  on  the  relative  time 
spent  on  the  two  dynamics.  Further  study  is  necessary  to: 

•  Optimize  the  decoupled  approach,  in  particular  with  respect  to  the  selection  of  the  tolerances  and  the 
step  sizes  internally. 

•  Test  the  decoupled  approach  using  more  complete  models  for  the  orbital  and  attitude  dynamics.  This 
will  also  allow  having  a  realistic  estimation  of  the  time  spent  on  the  evaluation  of  each  one  of  the 
dynamics. 

•  Include  three  dynamics  in  the  decoupled  approach  (orbit,  attitude  and  deformation),  or  alternatively 
embed  the  deformation  dynamics  in  the  attitude  dynamics. 
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8.  Appendix  A  -  ABAQUS  modelling  for  validation 


In  this  section,  we  use  finite  element  method  (FEM)  to  build  a  model  simulating  the  membrane  structure  at  GEO 
and  analysis  its  dynamic  behaviors  under  the  load  of  SRP.  As  for  material  properties,  we  used  the  most  common 
membrane-type  space  debris  the  separated  multi-layer  insulation  (MLI).  We  used  a  commercial  FEM  software 
ABAQUS  to  do  most  of  the  analysis  job  and  get  the  result,  some  of  which  were  used  as  a  validation  of  the  MATLAB 
of  the  deformable  debris  object. 

8.1.  Basic  ideas  of  FEM 

8.1.1.  Classification  and  application 

The  FEM  is  a  very  powerful  tool  for  solving  numerous  engineering  problems  and  widely  used  in  mechanics 
analysis,  thermal  analysis  and  many  other  branches  of  industry.  In  this  report  we  focus  on  mechanics  analysis. 

Considering  the  fact  that  load  applied  on  the  object  of  study  vary,  people  sort  the  problem  into  static  problem  and 
dynamic  problem,  then  use  different  governing  equations  to  describe  the  problem.  In  static  problem,  a  structure  is 
under  static  load  and  usually  reaches  equilibrium  condition.  Typical  static  problem  includes  structure  distortion 
analysis,  certain  structure  component  stress  and  strain  condition  analysis,  etc.  [51].  In  dynamic  problem,  however, 
the  direction,  value  and  application  point  of  the  load  can  all  be  functions  of  time,  and  the  FEM  is  used  to  calculate 
the  forced  response  of  the  body.  Typical  dynamic  problem  includes  forced  frequency  analysis  of  aircraft  and 
spacecraft  components,  analysis  of  impact  problem  and  time  dependent  analysis  of  structure  dynamic,  etc.  [51]. 

8.1.2.  General  theory 

The  basic  principles  underlying  the  finite  element  method  are  not  complex.  When  analyze  the  distribution  of  an 
unknown  variable  on  an  object,  taking  displacement  or  stress  in  elasticity  problem  as  an  example,  we  follow  the 
steps  below. 

The  first  step  is  discretization  of  the  problem.  The  real  object  with  infinite  degree  of  freedom  is  divided  into  a 
series  of  regular  shaped  and  manageable  pieces  called  “elements”,  and  they  are  interconnected  to  each  other  at  joints 
called  “nodes”.  In  FEM  every  element  represents  a  discrete  portion  of  the  real  object,  and  our  interested  variable, 
such  as  displacement,  is  calculated  only  at  the  nodes  of  the  element.  Then,  at  any  other  point  in  the  element,  the 
displacements  are  obtained  by  interpolating  from  the  nodal  displacements  in  a  certain  manner,  such  as  linear  or 
quadratic  polynomial  and  trigonometric  function,  which  is  usually  determined  by  the  element  type  we  choose  and  the 
number  of  nodes  used  in  the  element.  Therefore,  once  we  get  the  value  of  the  variable  at  every  node  of  every 
element,  the  variable  distribution  through  every  element,  and  further  the  whole  body,  is  adequately  approximated. 

All  structures  in  real  world  are  three-dimensional,  however,  considering  memory  usage,  computing  time  and  the 
property  of  the  problem  itself,  proper  simplification  is  essential  and  reasonable.  For  example,  it  is  economical  to 
develop  a  kind  of  element  that  can  carry  only  tensile  or  compressive  loads  without  resistance  to  bending  and  have 
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only  translational  degrees  of  freedom  at  each  node,  to  discretize  a  long,  axisymmetric  part  under  only  axial  load,  like 
a  part  of  the  scaffold.  Hence,  different  types  or  families  of  element  were  developed  for  different  simulation  and  result 
demand,  and  the  truss  element  is  exactly  the  element  type  stated  in  the  example  above.  Fig.  70,  from  [52]. 


Continuum 
(solid)  elements 


r\ 


Beam 

elements 

JAA/V 

* - E— 


Springs  and  dash  pots 


Truss 

elements 


Membrane  Infinite 

elements  elements 

Fig.  70.  Most  element  families  currently  available  in  commercial  FEM  software  for  mechanic  analysis  (from 
[52]) 


With  different  element  families  to  choose  and  adjustable  nodes  to  put  on  a  certain  element,  the  discretization  job 
becomes  more  flexible.  In  mechanic  simulation,  the  fundamental  variables  are  the  translational  degrees  of  freedom 
(DOF)  in  three  directions  and  rotational  DOF  in  three  directions.  In  order  to  get  the  required  variables,  such  as  stress 
strain  or  deformation  of  the  object,  they  should  be  calculated  step  by  step,  node  by  node  during  an  analysis  producer. 
Some  element  families  developed  for  certain  kinds  of  problems,  however,  eliminate  some  unnecessary  DOF  to 
simplify  the  model  and  accelerate  calculation.  What  if  the  element  family  is  suitable  but  a  more  accurate  distribution 
of  the  variable  is  desired?  A  more-node-type  element  is  developed  to  increase  the  interpolation  order  without 
changing  other  properties  of  the  same  element  family,  as  illustrated  in  the  examples  in  Fig.  71  (from  [52]) 


(a)  Lin  ear  etement 
(B-node  brick,  C3D8) 

Fig.  71.  Finite  element  types  (from  [52]). 


(b)  Quadratic  element 
(30-node  brick,  C3D20) 


•  Elements  that  have  nodes  only  at  their  corners,  such  as  the  8 -node  brick  shown  in  Figure  2(a),  use  linear 
interpolation  in  each  direction  and  are  often  called  linear  elements  or  first-order  elements. 

•  Elements  with  mid-side  nodes,  such  as  the  20-node  brick  shown  in  Fig.  71  (b),  use  quadratic 
interpolation  and  are  often  called  quadratic  elements  or  second-order  elements. 
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The  next  step  after  the  discretization  is  to  set  up  the  basic  governing  equation,  usually  called  finite  element 
equation.  Theoretically,  the  finite  element  equation  should  contains  all  the  information  of  the  system  that  we  know 
and  want  to  know,  for  instance  applied  load,  boundary  condition  and  displacements.  Finite  element  equations  differ 
in  formulation  for  different  problems,  however,  the  principle  to  derive  it  is  the  same,  to  set  up  the  equations  on  every 
single  element,  and  then  assemble  them  to  a  whole. 

Our  problem  belongs  to  the  elasticity  problem.  In  elasticity  analysis,  with  the  principle  of  superposition  and 
several  basic  type  of  load,  namely  concentrated  force,  distributed  force  or  pressure,  body  force  (inertia  or 
gravitational)  and  initial  strain  gO,  any  complex  force  condition  can  be  imitated  and  then  substituted  into  the  finite 
element  equations  to  yield  the  nodal  displacements.  Then  the  strain,  stress  and  other  variables  in  elasticity  analysis  is 
easily  to  get  from  the  nodal  displacements. 

The  derivation  of  the  general  finite  element  equation  for  elasticity  problem  is  given  by  the  method  of 
minimization  of  the  potential  energy. 

The  potential  energy  II  equals  the  strain  energy  A  in  the  body  less  the  work  done  by  the  external  loads  acting  on 
the  body: 

n  =  a  -  w 

In  a  differential  element  of  volume  d  V  ,  the  strain  energy  d  A  is  found  from 

dA  =  —  {s}T  {a}  -  —  {£0}T  {a}  (6) 

2  2 


Where,  {s}  is  a  vector  of  the  total  strain;  {^0}  is  a  vector  of  initial  strain;  and{cr}  is  a  vector  of  the  stress 
components. 

The  total  strain  energy  for  a  finite  volume  is  therefore  calculated  by  integrating  (6)  though  the  volume.  Hence 

A  =  f  ~({s}T{G}-{E0}T{G})dV  (7) 

J  v  2 


According  to  basic  elastic  theory,  the  strain  { s }  and  stress  { cr }  have  relation  as  below 
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(8) 


Where  [D]  is  elastic  coefficient  matrix,  determined  by  material  properties,  like  Young’s  modulus  E  and  Poisson 
Ratio  v  . 

Therefore,  if  we  substitute  (8)  into  (7)  we  can  get 

A  -  ^({£}r  [£>]{£}-  2{8}r[D]{80}  +  {EQ}T[D]{E0})dV 


Eql 


98 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Dynamics  of  flexible  MLI-type  debris  for  accurate  orbit  prediction 


«B  B3  University  School  of 
^  of Glasgow  Engineering 


Another  basic  principle  from  elastic  theory  is  that  the  strains  are  partial  differentiations  of 

T 

in  three  directions. 


displacements { u }  =  [ux  u2  u3 
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As  mentioned  before,  in  FEM  the  displacements  throughout  each  element  are  interpolated  from  the  nodal 
displacements  in  certain  manners. 

{u}  =  [NHU  } 

Where  { u  }  will  represent  the  displacements  of  any  point  on  an  element,  [  N  ]  is  the  shape  function  matrix  determined 
by  the  interpolation  function  and  {U  }  is  a  row  vector  of  nodal  displacements. 

Then  the  strains  can  also  equated  to  the  nodal  displacement  by 

{8}  =  [B]{U} 

Eq  2 

Compare  Eq  1  and  Eq  2,  we  can  find  that  [B]  is  a  matrix  derived  by  the  proper  differentiation  of  the  shape  function 
matrix  [  N  ]  . 

Substitution  of  Eq  2  into  Eq  1  gets  the  strain  energy  for  a  single  element  A (e)  ,  namely 


A  (e) 


-c 


({U  U)}T  [B^e)f  [D(e)][B(e)]{U  (e)}  - 


2{£/  (e  }  [ B(e) ]  [D^KSol  +  fBo  }  [D{e)]{z(0e)})  dV 


Eq  3 

In  order  to  get  the  potential  energy  of  a  single  element,  the  work  done  on  the  single  element  by  external  load  is  then 
needed.  As  mentioned  before,  no  matter  how  complex  a  real  load  may  be,  its  work  can  be  simplified  into  three  basic 
parts  that  are  the  work  done  by  concentrated  nodal  force,  by  pressure  and  by  body  force.  All  those  work  is  simply  the 
magnitude  of  the  forces  multiplying  the  distance  moved  on  respective  direction. 

Therefore,  if  the  vector  of  the  concentrated  nodal  loads  is  {  f  ie) }  ,  then  the  work  done  should  be 


=  {[/  U)j7-{F(e)j 


Eq  4 

T 

Assume |  p(xe)  p(ye)  p[e)}  is  the  distributed  pressure  loads  parallel  to  the  coordinate  axes  and  considered  positive 
when  acting  in  the  positive  coordinate  direction,  then  the  work  done  by  pressure  is 

r^ri 

T  |  | 

=  {  {iy (e)} [iv <e)]  {p\e)\ds 


(e) 

[Pz  J 


Eq  5 


Assume  the  body  force  isjx(e)  Y{e)  Z  (e)  j  per  unit  volume,  then  the  work  done  by  body  force  i 


is 
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fx<e,l 


W[Be]  =r  {U(e)}T[N(e)]T  {Ym  \dV 

L(«)  l 

lz  J 

Eq  6 

Therefor  the  potential  energy  of  a  single  element  is  given  by  Error!  Reference  source  not  found,  and  Eq  3  to  Eq  6, 
namely 

n <c>  =  a  <e)  - 1 yc(e)  -  w -  wBe) 

Eq  7 

If  the  real  object  is  discretized  into  N  elements,  then  the  potential  energy  of  the  whole  system  is  simply  the  sum  of 
every  single  element 

N 

n  =  £  n  M 

e  =  \ 

Eq  8 

For  the  object  under  load  to  be  in  equilibrium,  its  potential  energy  must  reach  a  minimum.  Hence  for  all  degrees  of 
freedom,  namely  displacements  for  every  node  on  all  three  possible  directions,  must  obey 

an  an  an 

du ,  du _  du 

12  n 


or 


an 


=  0 


Eq  9 

Where  we  assume  n  is  the  total  degrees  of  freedom,  {u  s }  is  the  vector  consists  of  all  elemental  displacement  vectors. 
Calculate  out  Eq  9  we  can  actually  get  n  equations  with  n  unknowns 

^  j-j  N 

- =  Z  [  f  ,  ,  ([S<e)]r  [Z3  “,)  ][s“'1])  dV  ]{t/si 

5{t/  ^  J''<> 


fz'-M 

-2  ([JBt',]r[o<')]{s^,i)  dV  +  jr1*'  \dv 

\  ry  (.)  I 

lz  J 

\pM  1 

T  I  I 

+  <‘°}  [AT00]  ^  p\e)  \dS  ]-{PJ  =  0 


To  be  simplify  Eq  10  may  be  written  as 


Eq  10 


[K]{US}  =  {FS] 


Eq  11 


Where 

N  N 

IK]=  X  [Ji/(„)([5<e>]r[£»<e)][S<e)])  dV]  = £  [*<e)] 

e=\  e=l 

Eq  12 

[K  ]  is  global  stiffness  matrix,  [K  {e)  ]  is  the  individual  stiffness  matrix  and  the  global  force  vector  is 
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fxU)] 

N  |  | 

{Fs}  =  Z  [  Jv,.,  dV  +  J  U<e)  \dV 

-1  I  _(.)  I 

lZ  J 

r  I  I 

+  f  {[/<e>}[iV<e)]  ]  +  {Ps}  =  Z 

I  I  <•=> 

L p*  J 

Eq  13 

Usually  the  individual  element  stiffness  matrix  and  force  vector  for  each  element  is  solved  first,  and  then  they  are 
assembled  into  the  global  matrix  and  vector.  Certainly  the  symbol  “X”  does  not  mean  what  it  usually  means.  The 
individual  stiffness  matrixes  are  obtained  in  the  local  coordinate  of  every  element,  to  assemble  into  the  global  matrix 
they  should  first  be  converted  into  global  coordinate  and  then  added  to  the  proper  locations  in  the  global  matrix 
according  to  its  nodal  association. 

Now  we  have  derived  Eq  1 1  as  the  basic  equation  for  finite  element  problem.  It  is  the  governing  equation  for  elastic 
static  analysis,  however  for  dynamic  analysis,  which  is  exactly  our  situation,  it  is  not  comprehensive  enough.  There 
should  be  some  other  terms  in  the  equation  that  can  represent  the  features  of  dynamic  problem,  namely  inertia  and 
damping.  So  the  basic  equation  of  motion  for  the  object  becomes 


[M  ]{t/5}+[C]{U5}+[^]{C/5}  =  {Fs(t)} 


In  which,  [  M  ]  is  mass  matrix  and  [  C  ]  is  damping  matrix,  which 


Eq  14 

are  determined  by  the  material  properties, 


and { Fs  (t) }  means  the  load  may  vary  as  time  goes  by. 

After  we  have  get  the  finite  element  equation  for  the  problem  the  final  step  is  to  solve  the  equation,  there  are  mainly 
two  classes  the  direct  methods  and  the  iterative  methods,  both  of  which  have  been  studied  a  lot  in  mathematics  and 
numerical  analysis,  and  most  commercial  FEM  packages  have  integrated  the  latest  achievements.  We  will  not  focus 
on  the  solution  methods  in  this  report. 


8.1.3.  A  brief  introduction  to  ABAQUS 


ABAQUS,  based  on  FEM,  is  a  suite  of  engineering  simulation  programs,  including:  ABAQUS/Standard,  a 
general-purpose  finite  element  program;  ABAQUS/Explicit,  an  explicit  dynamics  finite  element  program; 
ABAQUS/CAE,  an  interactive  environment  used  to  create  finite  element  models,  submit  analyses,  monitor  and 
diagnose  jobs,  and  evaluate  results;  and  ABAQUS/Viewer,  a  subset  of  ABAQUS/CAE  that  contains  only  the  post¬ 
processing  capabilities  of  the  Visualization  module.  It  is  a  general-purpose  simulation  tool  that  can  be  used  in  diverse 
areas,  such  as  stress/strain  analysis,  heat  transfer  analysis,  mass  diffusion  analysis,  piezoelectric  analysis, 
electromagnetic  analysis,  and  fluid  dynamics  [6].  In  ABAQUS/CAE  there  is  an  extensive  elements  library  to  model 
complex  geometry  and  an  extensive  material  library  to  simulate  the  behavior  of  most  typical  engineering  materials. 

Using  ABAQUS  to  analyze  a  problem  should  follow  three  successive  steps,  preprocessing,  simulation,  and  post¬ 
processing,  shown  in  Fig.  72  three  steps  to  analyze  problem  by  ABAQUS  [6]: 
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Fig.  72  three  steps  to  analyze  problem  by  ABAQUS 

A  basic  concept  in  ABAQUS  is  the  division  of  the  problem  history  into  “step”  that  is  any  phase  of  the  history. 
Taking  the  simplest  form  in  ABAQU S/Standard  static  analysis  as  an  example,  a  “step”  is  a  load  change  from  one 
magnitude  to  another.  In  each  “step”  the  user  chooses  a  “procedure”,  defining  the  type  of  analysis  to  be  performed 
during  the  step:  static  analysis,  dynamic  stress  analysis,  eigenvalue  buckling,  etc.  The  procedure  choice  can  be 
changed  from  step  to  step  in  any  meaningful  way.  Since  the  state  of  the  model  is  updated  throughout  all  analysis 
steps,  the  effects  of  previous  steps  are  always  included  in  the  response  in  each  new  step. 

ABAQUS/Standard  provides  both  linear  and  nonlinear  response  options.  Linear  analysis  is  considered  as  linear 
perturbation  analysis.  This  linear  perturbation  approach  allows  general  application  of  linear  analysis  techniques  in 
cases  where  the  linear  response  depends  on  preloading  or  the  nonlinear  response  history  of  the  model  [6] . 

In  nonlinear  analysis,  which  is  our  case,  the  goal  is  to  obtain  a  convergent  solution  by  a  proper  method  at  the 
minimum  cost.  So  the  solution  will  be  developed  by  a  series  of  “small”  increments.  For  a  successful  FEM  package, 
two  aspects  must  be  focused  on:  first,  the  method  solving  discrete  equilibrium  equation  at  each  increment;  second, 
the  technique  to  setting  the  proper  increment  size.  For  solving  the  nonlinear  equilibrium  equations  at  each  increment, 
many  numerical  methods,  usually  based  on  Newton  method,  are  available  for  ABAQUS/Standard.  Again  taking  our 
problem  as  an  example,  Modified  Newton's  method  is  just  suitable.  Other  methods  like  Full  Newton  method, 
Approximate  Newton  method  and  Quasi-Newton  method,  etc.  all  have  their  suitable  cases. 

While  setting  the  increment  size,  certainly  ABAQUS  will  allow  the  direct  user  control  of  increment  size  where 
the  user  specifies  the  incrementation  scheme.  The  automatic  control  of  increment  size,  however,  should  manifest  the 
value  of  ABAQUS  more.  This  approach,  in  which  user  defines  the  step  and  specifies  certain  tolerances  or  error 
measures  and  then  let  ABAQUS  automatically  selects  the  increments  as  it  develops  the  response  in  the  step,  is 
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usually  more  efficient,  since  it  is  hard  for  user  to  predict  the  response  ahead  of  time.  Automatic  control  is  particularly 
valuable  in  cases  where  the  time  or  load  increment  varies  widely  through  the  step,  such  as  creep,  heat  transfer,  and 
consolidation  [6]. 

Those  are  just  brief  introduction  of  ABAQUS,  more  detailed  discussion  on  solving  method  and  theory  shall  refer 
to  ABAQUS  user’s  Manual  or  other  reference  literature. 

8.2.  Modeling  and  Result 

In  this  project,  our  work  was  built  in  ABAQUS/CAE  and  use  ABAQUS/Standard  to  carry  out  the 
stress/displacement  analysis.  And  the  properties  we  used  to  model  the  space  membrane  structure,  as  we  have 
mentioned  in  the  first  part  of  the  report,  are  from  a  sheet  of  separated  MLI  and  shown  in  the  table  below. 


Young’s  modulus  :  8.81  *  109  N/m2 


Poisson’s  Ratio  :  0.38 


Side-length  :  lm 


Thickness  :  2.5xl0'6m 


Density  :  1.39xl03  kg/m3 


Solar  pressure  :  8.01  lxlO'9  pa 

(When  the  face  normal  directed  at  the  sun) 


In  our  simulation  there  is  a  thin  membrane  under  relatively  small  load  changing  its  geometry  slowly,  and  we  are 
interested  in  the  distribution  of  displacement,  stress  and  strain.  So  based  on  these  characteristics  of  our  problem,  the 
element  type  we  chosen  is  linear  (4-node)  quadrilateral  curved  shell.  It  is  a  general-purpose  element  that  can  provide 
robust  and  accurate  simulation  for  thin  and  or  shell  under  various  loads.  When  using  this  element,  the  membrane  is 
assumed  having  no  transverse  shear  locking,  and  unconstrained  hourglass  modes.  This  element  considers  finite 
membrane  strains  and  the  membrane  kinematics  are  based  on  an  assumed -strain  formulation  that  provides  accurate 
solutions  for  in-plane  bending  behavior. 

There  are  three  typical  original  membrane  geometries  designed  to  analyze  its  response  under  load  in  space 
(mainly  solar  pressure). 

8.2.1.  Case  1,  Flat  plate  under  normal  pressure 

Assume  the  face  normal  directed  at  the  sun  and  the  solar  pressure  applied  for  1  second,  shown  in  Fig.  73,  the  pink 
arrow  represents  the  force  vector. 


103 

Distribution  A:  Approved  for  public  release;  distribution  is  unlimited. 


Dynamics  of  flexible  MLI-type  debris  for  accurate  orbit  prediction 


|g  University  School  of 
of Glasgow  Engineering 


And  here  is  the  simulation  result  of  displacement  from  ABAQUS 


As  we  can  see  from  Fig.  74  the  whole  body  is  in  the  same  color,  which  means  the  displacement  of  every  element 
is  the  same,  as  it  marked  out  in  white,  max=min=1.152e'6  (m).  As  for  the  stress  and  strain,  it  is  no  surprise  to  find 
they  are  zero  all  around  the  plate. 
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Actually,  for  this  simplest  situation,  the  solar  pressure  is  added  evenly  to  the  surface,  it  is  not  hard  to  get  the  same 


conclusion  based  on  Classical  Kinematics  Theory.  For  all  elements,  they  have  the  same  mechanics  properties  and  the 


load  situations,  so  after  the  period  of  time  they  should  have  the  same  movements.  And  the  displacements  of  every 
element  can  be  calculated  out  equally  by  Kinematics  theory. 

From  Newton's  laws  of  motion  the  acceleration  equals: 

F 

a  =  — 

m 

where:  F=8.01  lxlO'9xlx  1=8.01  lxlO'9  N,  is  the  total  load  applied  on  the  body 
m=2.5xl0'6xlxlxl.39xl03=3.475xl0'3  kg,  is  the  total  mass  of  the  body. 

So,  a=2.305xl0~6  m/s2 

Then,  from  the  basic  Kinematics  equation: 


the  displacement  d=  1.1 5 2x1  O'6  m. 

The  results  of  the  same  problem  from  two  different  methods  are  same,  which  to  a  certain  extent  showed  the 
correctness  of  this  method. 

8.2.2.  Case  2,  Right-angle  membrane  under  constant  solar  pressure 

In  this  case  we  assume  the  original  geometry  for  this  membrane  is  not  plate.  As  we  can  see  in  Fig.  75  it  forms  a 
right-angle  at  the  middle  point  of  edges.  The  solar  pressure,  added  to  the  face-YOZ,  will  last  for  60  seconds,  and  to 
make  the  simulation  more  factual  its  direction  is  defined  with  respect  to  the  original  body,  which  means  the  load 
always 


Fig.  75  force  conditon  for  case  2 


Here  are  the  simulation  results  for  this  case.  Fig.  76  shows  displacement  distribution  of  every  node  after  60 


seconds.  As  it  marked  out,  the  maximum  and  minimum  displacement  located  respectively  at  the  upper  edge  and 
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lower  edge,  max=3. 822x1  O'3,  min=1.043xl0"3.  Fig.  77  and  Fig.  78  shows,  respectively,  the  distribution  of  stress 
and  logarithmic  strain  at  60  second,  also  marking  out  the  max  and  min  values,  and  still  their  contours  pictures  are 
same. 


Fig.  76.  displacement  distribution  for  case  2 


Fig.  77  stress  distribution  for  case  2. 
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Fig.  78  strain  distribution  for  case  2 


In  real  situation  the  membrane  is  under  a  very  complex  force  condition  and  the  force  applied  on  it  may  change  as 
time  goes  by.  So  in  the  following  case  we  simulate  the  situation  that  the  membrane,  in  the  same  origin  geometry  of 
last  one,  was  destabilized  by  a  relatively  large  force  in  the  first  five  seconds  of  the  total  duration  of  solar  pressure. 


8.2.3.  Case  3,  Right-angle  membrane  under  changeable  load 

In  this  case,  we  assume  in  the  first  five  seconds,  beside  the  solar  pressure,  there  is  another  pressure  applied  on  the 
membrane.  As  shown  in  Fig.  79,  the  new  arrows  normal  to  the  vertical  half-face  represent  the  new  pressure  and  its 
magnitude  is  10'5  pa. 


Fig. 


79  force  conditon  for  case  3 
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Fig.  80  shows  the  displacement  distribution  under  the  new  disturbing  load,  as  we  can  see,  the  transparent  body  is 
the  original  body  and  the  entity  is  how  it  deformed  after  60  seconds.  Obviously,  this  distinguishable  deformation  will 
have  a  complex  influence  on  its  motion  characteristic.  Fig.  81  is  the  stress  contours  picture  for  this  case,  and  still  the 
strain  has  the  same  contours  picture  as  stress,  which  is  Fig.  82. 


Fig.  80  displacement  distribution  for  case  3 


Fig.  81  stress  distribution  for  case  3 
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Fig.  82  strain  distribution  for  case  3 


8.2.4.  Case  4 ,  Curve  membrane  under  constant  solar  pressure 

In  this  case  the  original  geometry  is  semi-cylinder,  as  we  can  see  in  Fig.  83.  The  solar  pressure,  lasting  for  60 
seconds,  is  applied  to  the  membrane  from  negative  Y  direction,  defined  with  respect  to  the  original  body. 
Considering  the  essence  of  solar  pressure  is  particles  collision,  it  is  no  surprise  having  the  load  distribution  shown  in 
Fig.  83. 


Fig.  83  force  conditon  for  case  4 

The  following  three  Figures  show  the  simulation  results  for  this  case.  Fig.  84  shows  the  displacement  distribution 
after  60  seconds.  As  it  marked  out,  the  maximum  and  minimum  displacement  is  max=2.907xl0'3,  min=2.496xl0'3, 
respectively.  Fig.  85  and  Fig.  86  shows  the  distribution  of  stress  and  strain  at  60  second,  respectively,  also  marking 
out  the  max  and  min  value. 
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Fig.  84  displacement  distribution  for  case  4 
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Fig.  86  strain  distribution  for  case  4 

8.2.5.  Case  5,  Curve  membrane  under  changeable  load 

Like  Case  3,  we  assume  in  the  first  five  seconds,  beside  the  solar  pressure,  there  is  another  face  load,  having  the 
same  characteristic  of  solar  pressure,  applied  on  the  membrane  from  negative  X  direction,  as  shown  in  Fig.  87,  and 
its  magnitude  is  10'6  pa. 


Fig.  87  force  conditon  for  case  5 

Similarly,  Fig.  88  shows  the  displacement  distribution  in  60  seconds  for  this  case,  max=0.0502  and  min=0.0102. 
Fig.  89.  and  Fig.  90  is  the  contours  picture  for  stress  as  well  as  strain.  Like  Case  3,  with  a  relatively  large  initial 
disturbing  force,  the  membrane  will  perform  a  more  obvious  deformation,  which  in  result  will  lead  to  larger  and 
more  complex  displacement. 
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Fig.  88.  displacement  distribution  for  case  5 
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Fig.  90  strain  distribution  for  case  5 


8.3.  Conclusion 

In  this  report  we  use  the  properties  of  MLI  to  build  a  model  for  a  simple  square,  thin -membrane  structure.  Then  three 
origin  geometries  and  five  force  conditions  was  set  as  a  simulation  of  real  cases.  For  every  case,  we  use  ABAQUS  to 
simulate  the  motion  the  deformation  and  display  the  result  in  contours  pictures. 

From  those  results  we  can  find  in  space  even  very  small  force  lasts  for  very  short  period  will  have  a  nonnegligible 
influence  on  the  thin  flexible  membranes,  especially  considering  their  deformation.  The  very  easily  happened 
deformation  will  change  their  load  condition  and  in  turn  change  the  motion  law  of  those  on  orbit  membranes. 
Another  interesting  finding  is  the  contours  pictures  for  stress  and  strain  is  always  the  same  which  means  the 
deformed  structure  will  always  have  the  same  distribution  on  stress  and  strain. 

Actually,  what  we  have  done  now  is  just  a  starter  of  the  study.  In  order  to  build  up  an  accurate  model  to  describe  the 
deformation  as  well  as  attitude  and  orbit  changing  of  on  orbit  membrane,  there  is  a  lot  more  need  to  improve  and 
prefect.  As  we  know  it  is  a  long  haul,  usually  measured  by  hours  or  days,  that  complex  loads  in  space  influence  the 
on  orbit  membrane.  In  our  analysis,  however,  the  duration  of  the  load  is  quite  short  and  the  load  condition  we 
considered  is  also  simplified,  mainly  solar  pressure  and  a  simple  disturbing  force  in  two  cases.  Therefore,  in  next 
step,  it  is  necessary  to  lengthen  the  load  duration  and  build  a  more  accurate  load  condition  model  based  on  the 
location  of  the  membrane  on  the  orbit.  Besides,  more  complicated  structure  and  origin  geometries  should  be 
considered. 
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